NASA Technical Memorandum 100163 

4» ' 


Numerical Simulation of Subsonic 
and Transonic Propeller Flow 


(KASA-TM-100U3) NUMERICAL SIMULATION OF 
SUBSONIC AND TfiANSCNIC PBOFBIIEB FLOW Ph.D. 
liesis (NASA) 191 p CSCL 01A 


G3/02 


Aaron Snyder 
Lewis Research Center 
Cleveland, Ohio 


April 1988 


N88-2 C2€2 

Uuclas 

0133311 




NUMERICAL SIMULATION OF SUBSONIC AND TRANSONIC PROPELLER FLOW 

Aaron Snyder 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 


SUMMARY 


The numerical simulation of three-dimensional transonic flow about 
a system of propeller blades Is Investigated. In particular, It Is 
shown that the use of helical coordinates significantly simplifies the 
form of the governing equation when the propeller system Is assumed to 
be surrounded by an Irrotatlonal flow field of an Invlscld fluid. The 
unsteady small disturbance equation, valid for lightly loaded blades and 
expressed In helical coordinates, Is derived from the general blade- 
fixed potential equation, given for an arbitrary coordinate system. The 
use of a coordinate system which Inherently adapts to the mean flow 
results In a disturbance equation requiring relatively few terms to 
accurately model the physics of the flow. Furthermore, the helical 
coordinate system presented here Is novel In that It Is periodic In the 
circumferential direction while, simultaneously, maintaining orthogonal 
properties at the mean blade locations. The periodic characteristic 
allows a complete cascade of blades to be treated, and the orthogonality 
property affords straightforward treatment of blade boundary conditions. 
An ADI numerical scheme is used to compute the solution to the steady 
flow as an asymptotic limit of an unsteady flow. As an example of the 
method, solutions are presented for subsonic and transonic flow about a 
5 percent thick bicircular arc blade of an eight bladed cascade. Both 
high and low advance ratio cases are computed, and include a lifting as 
well as nonlifting cases. The nonlifting solutions obtained are 
compared to solutions from an Euler code. 
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I. INTRODUCTION 


The rap 1 c J rise of computer capability has greatly influenced the 
methodology for attacking complex problems in engineering and science. 
The trend toward the numerical formulation and solution of problems has 
grown and will undoubtedly expand in the future. In the areas of fluid 
mechanics and heat transfer, this approach is known as computational 
fluid dynamics (CFD) . 

The application of computational fluid dynamics methods toward the 
simulation of the flow surrounding high-speed propellers is the subject 
of this work. 

Many examples of both steady and unsteady flows of practical 
interest are well approximated by various forms of the governing 
equations of fluid mechanics. The recent marked advance in computer 
capabilities and the rapid growth of efficient and accurate numerical 
methods have allowed a number of difficult flow problems to be examined. 

In the past, analytical and experimental methods were the two main 
avenues for obtaining aerodynamic solutions. While analytical methods 
are still being enhanced and often give quick solutions, they suffer the 
drawback of being limited to simple configurations. On the other hand, 
wind-tunnel experiments can handle complicated configurations, but are 
limited to certain flow ranges and require expensive models and 
extensive test times. Nevertheless, the wind tunnel has served as the 
primary tool for the development of aerodynamic configurations. Shapes 
can be tested and modified in light of pressure and force measurements. 
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Analytical methods are generally limited to simpler isolated components 
and partially owe their usefulness to the fact that efficient 
aerodynamics can only be achieved with well ordered flows. 

Consequently, simplified mathematical models are often adequate 
replacements for the full set of governing equations. This was utilized 
with great success in many important applications. For example, often 
purely inviscid flow assumptions for irrotational flow, in which the 
velocity c?n be formulated in terms of a gradient of a scalar, have 
proven sufficient. However, both the limitations of the analytical 
techniques as well as the costs of the wind tunnel studies are all too 
apparent as the features of the flow fields become more complex. This 
is especially true for the prediction of transonic flows. Thus, a need 
for alternative solution methods was created. 

The growing area of computational fluid mechanics offers an 
attractive complementary alternative to analytical methods or 
experimentation by providing details of the flow physics as well as data 
beyond the experimental range while affording configuration 
optimization. The use of computational methods to provide solutions 
about rotating propellers is one illustration of the ability of these 
techniques to solve a complicated set of nonlinear partial differential 
equations and associated boundary conditions. 

Coinciding with the arrival of increasingly powerful CFD methods, 
was the emergence of renewed interest in propellers which offer large 
performance improvements for aircraft that cruise ir. the high subsonic 
speed range. The reason for this is that for lifting surfaces, to a 
first approximation, the efficiency is proportional to the Mach number 
of the approaching flow times the lift to drag ratio. The lift to drag 
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ratio rises monotoni cal ly with Mach number until a sharp drop-off occurs 
at high subsonic free-stream Mach numbers. Thus, it pays to increase 
the speed as much as possible until losses associated with increased 
drag become prohibitive. 

The very high propulsive efficiencies of the advanced free air 
propeller coup’ed with the need to reduce fuel costs has led to a 
rebirth in propeller research programs. These programs have introduced 
changes that have altered advanced high-speed turboprops considerably 
from their low-speed turboprop counterparts. The most notable changes 
are the use of eight to ten thin, highly swept, smal 1 -diameter blades 
instead of two to four longer, thicker, and unswept conventional blades. 
These changes allow propellers to operate efficiently in the transonic 
flow regime. However, they require sophisticated computational methods 
in order to predict and calculate their aerodynamic characteristics. 

This is mainly due to the increased three-dimensionality of the flow 
field encompassing these low aspect ratio blades and to the rise of 
mathematical nonlinearity in the flow model which is associated with 
their transonic tip speeds. Also, their thin structural design makes 
the blades prone to flutter. This is a source of unsteadiness which 
must be considered in addition to the inherent unsteadiness of the 
transonic regime. 

Improvements in present techniques are needed to provide accurate 
solutions for unsteady three-dimensional transonic flows about 
propellers. Tne use of numerical methods offers the potential for 
solving complex equations in intricate geometries such as those which 
exist in the case of simulating flow about rotating propeller blades. 
However, as with any method, limits exist on the capabilities. These 
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limits involve the speed and memory of the computers and the 
availability of efficient algorithms as well as considerations of 
robustness and affordability. For a given problem, judgment dictates 
the trade-off between the suitability of a particular model and the 
feasibility of the calculation. This is very much the same dilemma as 
faced when using a theoretical approach. 

It is generally accepted that the full Navier-Stokes equations, 
which satisfy conservation of mass, momentum and energy in a viscous 
heat conducting fluid, provide a complete description of most flows of 
interest to CFD. The unsteady compressible Navier-Stokes equations are 
a mixed set of hyperbolic-parabolic equations; for unsteady 
incompressiole flow they are a mixed set of elliptic-parabolic 
equations; for steady compressible flow, they are a set of hyperbolic- 
elliptic equat s ons . In the latter case, due to the effect of the 
different character of the numerical schemes required for hyperbolic and 
elliptic terms, the equations are often cast in the unsteady form and 
marched to a steady solution. At high Reynolds number, the unsteady 
compressible Navier-Stokes equations are difficult to solve due to the 
large difference in magnitude between the inertial and viscous forces. 
The large inertial forces of the hyperbolic terms can impose small time 
steps in order to meet accuracy or stability requirements. This time- 
step limitation generally severely retards the effects of viscosity 
associated with the parabolic terms of the equation. The consequence 
is that excessive computer time is needed to resolve the flow. 

The solution of the Navier-Stokes equations for general three- 
dimensional, time-dependent, turbulent flow is a distant goal. 

Candidates for the near future include certain averaged forms of the 
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Navier-Stokes equations. The one most likely to find application to 
aerodynamics Is the Reynolds-averaged form. Here the the turbulent 
terms are averaged with respect to time and the stresses are determined 
by empirical closure formulas. Other more advanced forms of modeling 
exist such as the large eddy simulation model. In this case the full 
time-dependent Navier-Stokes equation is used to calculate the 
development of large scale turbulent eddies. These eddies are 
anisotropic and highly structured. A process analogous to Reynolds 
averaging is used to model the subscale turbulence. Although 
interesting results have been reported, computational time for the 
solution of either the Reynolds averaged model or large eddy simulation 
model is exorbitant for general use. 

Current progress in aerodynamic design has relied mainly on the 
approach fostered by Prandtl, whereby viscous effects can be considered 
to be confined to a relatively thin region, called a boundary layer, 
adjacent to solid surfaces and their wakes; the entire remaining 
external flow can be regarded as being inviscid. This method 
effectively splits the flow into two distinct but coupled regions. The 
forms of the governing equations have quite different characteristics 
and demand different solution methods. This allows a considerable 
savings in computer time and storage since all the viscous terms are 
neglected in the inviscid region and only the viscous terms normal to 
the wall are retained in the boundary layer. In addition, the momentum 
equation is simplified in the viscous region and as a consequence the 
normal pressure gradient is neglected. Special boundary-layer type 
marching techniques are available which offer additional savings in 
computation time and storage. Much attention has been focused on 
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providing the correct representation of the vi scous-i nvi sci d coupling 
and the proper matching of their solutions at the interface of the 
regions. In many cases it is not acceptable to use two sets of 
equations. These cases include flows with strong coupling between the 
outer inviscid region and the boundary layer, such as found in a 
supersonic flow around a blunt body or in flows containing strong 
shocks. In these instances it would not be valid to neglect the normal 
pressure gradient when forming the normal momentum equation. As a 
consequence boundary layer techniques are often inadequate. 

Another set of equations which are approximations to the Navier- 
Stokes equations and which are valid for both the inviscid and viscous 
regions, are receiving considerable recent attention. These are known 
as the thin-layer Navi er-Stokes equations. They simplify the complete 
Navier-Stokes equations by neglecting the viscous terms containing 
derivatives parallel to a solid surface, as in the boundary layer 
formulation, but retain all the other terms in the momentum equation; 
thus, the pressure gradient in the direction normal to a surface is not 
neglected. The use of the thin-layer Navier-Stokes equations allows the 
robust calculation of separated and reverse flows and flows containing 
significant normal pressure gradients. In many instances the thin-layer 
equations can be treated by the same efficient boundary-layer type 
marching techniques that are employed in solving the boundary layer 
equations . 

Sets of equations, each of somewhat different form and similar to 
the thin layer Navier-Stokes equations, are obtained if only the 
streamwise viscous terms are omitted from the Navier-Stokes equations. 
These forms are called the parabolized Navier-Stokes equations, and the 
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differences in the forms depend on their derivation and whether terms 
such as the streamwise pressure gradient are retained in the streamwise 
momentum equation. The thin layer and the parabolized equations are 
equivalent for two-dimensional flows. However, efforts to develop 
efficient solution procedures are hampered by the fact that different 
forms generally require different strategies due to the nature of the 
equations or characteristics of the flow. 

Very useful inviscid forms of the governing equations exist. 
Currently, solutions of the Euler equations, obtained from the Navier- 
Stokes equations by neglecting all the viscous terms, are solved 
routinely for two-dimensional flows and even three-dimensional codes are 
emerging. Wider use has been made employing the full potential and 
small disturbance forms of the potential equation for inviscid, 
irrotational flows. While the Euler equations are valid for rotational 
flows, the potential formulations are not, and thus, the latter are not 
strictly valid when entropy producing discontinuities such as shocks 
exist. On the other hand, many instances occur where the departure from 
purely irrotational flow is slight, such as the flow around thin bodies 
or in the case of shocks which are sufficiently weak as to generate 
little entropy or vorticity. 

In general, considering the flow around an aerodynamic surface to 
be inviscid allows the lift to be calculated quite accurately. The 
drag, on the other hand, depends substantially on frictional stresses 
as well as the distribution of the normal pressure over the surface. 
Frictionless flow theory can be used for the determination of the drag 
component arising from the normal pressure components whereas an 
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appropriate skin friction coefficient must be used to obtain an estimate 
of the total drag. 

A number of numerical codes are currently being developed to 
calculate propeller flow. These solve either the full potential or 
Euler equations. Although these equations accurately model the flow, 
they are computationally more expensive and are more complicated to code 
than the small disturbance equation. 

A computer program was developed (Refs. 1 to 4) to simulate the 
flow over blade-tips of helicopter rotors. This code solves the small 
disturbance equation (Ref. 5) appropriate to a helicopter in forward 
flight. In particular, results were obtained that showed the code was 
able to track the development of a shock and its subsequent propagation 
upstream. The successful application of the small disturbance code in 
this manner indicated that the same approach could be applied to the 
propeller problem. 

The fast that turboprops operate at high subsonic cruise speeds 
means that at least the tip region of each blade will be embedded in a 
transonic helical flow. Shocks may likely exist ahead of, or on, the 
blades, depending on the flight speed and propeller design. Shocks of 
high strength indicate large wave drag. Efficient flight dictates that 
only weak shocks are acceptable, which implies little entropy 
generation. The potential formulation should provide a good 
approximation to the Euler equations as long as the shock strengths are 
less than that of a normal shock with an upstream Mach number less than 
about 1.3. Presently, turboprops that have helical Mach numbers 
slightly less than 1.0 at the hub to about 1.12 at the tip are being 
designed. It is reasonable to expect that the shock strength is 
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sufficiently weak near the blade tips, where three-dimensional relief 
effects occur. Hence, that the potential formulation is adequate. 
Furthermore, the use of a perturbation potential in place of a full 
potential significantly simplifies the formulation of the problem. 

Although many similarities exist between the equations developed 
for the helicopter rotor and the propeller problem, fundamental 
differences exist between the two. These differences demand that 
special consideration be made in accurately treating each type of flow. 
For example, in the case of the helicopter rotor, the flow field as 
observed In a coordinate system attached to the blade remains unsteady 
even if blade flutter is ignored. However, when the axis of rotation 
for a propeller coincides with its flight direction, a steady flow 
results with a blade attached coordinate system. This allows the steady 
problem to be examined separately from the flutter or gust problem. 

With the goal of using the numerical algorithm developed for the rotor 
code to the propeller problem, several concerns need to be addressed. 
Many of these concerns are naturally handled by using a suitable 
coordinate system. 

The key to simplifying the potential formulation of the propeller 
problem is the use of the helical free-stream direction as the primary 
direction. This provides an accurate primary flow upon which a 
perturbation can be superposed. The use of a helical flow, which 
inherently captures the flow curvature, results in a marked improvement 
over using an axial flow as the primary flow direction. The optimum 
flow to base a perturbation about is the exact flow; helical flow much 
more closely approximates the exact flow about a rotating and advancing 
propeller than does axial flow. 
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Since the helical flow direction captures the fundamental 
properties of propeller flow, relatively few terms are needed in the 
resulting disturbance equation to provide accurate modeling of the flow. 
In fact, no more terms are needed than when using Cartesian coordinates 
for the rotor problem. In addition, due to the similar nature of the 
two problems, the terms in the two sets of disturbance equations 
correspond identically except for their coefficients. Thus, the 
propeller equation is amenable to the same potential equation solver as 
developed for the helicopter. This is quite advantageous since much 
effort has been expended, both in making the potential solver efficient 
and in verifying its operation. 

While the helical flow direction is an ideal choice for the primary 
flow direction, other considerations must be taken into account in order 
to provide a satisfactory coordinate system. These considerations hinge 
on the desire to include cascade effects and the employment of small 
disturbance boundary conditions near the blade surfaces. Proper 
specification of the helical coordinates permits the straightforward 
treatment of annular cascade flows while maintaining orthogonal 
properties near the blade locations. The helical coordinates given in 
this thesis allow both of these requirements to be met. 

These helical coordinates are analytically defined in terms of 
basic propeller parameters through simple transformations to Cartesian 
coordinates. Thus, they automatically adapt to the cascade 
configuration for any number of blades and blade twist, which itself is 
a function of the propeller advance ratio. The helical coordinates 
enter the proolem in terms of the components of the metric tensor for 
the transformation. These components are readily evaluated in a 
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separate subroutine. The Jacobian of the transformation is particularly 
simple and provides a means of verifying whether a given transformation 
is valid. 

Obtaining accurate solutions for steady, not to mention unsteady, 
three-dimensional transonic flows, presents a difficult challenge even 
for the most advanced computational methods. The present work assumes 
a steady-state flow exists as observed in a reference frame attached to 
the propeller system, rotating with a constant angular velocity. This 
is considered a necessary step which must be completed before 
undertaking the unsteady problem. In addition, tne fluid will be 
regarded as being inviscid in an irrotational flow field and thus is 
expressible in terms of a potential formulation. The solution is sought 
in terms of a disturbance potential. This disturbance potential is the 
potential associated with the reduced velocity obtained by subtracting 
both the free-stream velocity and the rotational velocity, resulting 
from the transformation to a noninertial reference frame, from the total 
velocity. This solution should be valid for lightly loaded blades. 

Presentation of this thesis will begin with an introduction of some 
basic propeller concepts. In turn this will be followed by a historical 
review of some notable analytical and numerical endeavors related to 
propeller theory. This review is intended to provide the background and 
motivation for this work as well as form a basis for justification of 
the approach taken. In the subsequent chapter, the governing equation 
for the three-dimensional small perturbation potential is derived in 
detail from the general potential equation in a rotating reference 
frame. Here, the full potential equation is given in general tensor 
form applicable to any suitable coordinate system which need not be 



12 


orthogonal. Following this, the small disturbance form of the equation 
is derived in terms of coordinates which are scaled by the relevant 
transonic scaling parameters. A separate chapter is next provided for 
the boundary condition equations. A chapter dealing with the coordinate 
system construction and its arrangement to suit a cascade of propellers 
is then presanted. This is related to two appendixes on analytical 
helical coordinate transformations. Next, a chapter on the numerical 
procedure and solution technique is given. A generalized form of the 
Douglas-Gunn algorithm is introduced which allows for additional cross- 
derivative terms to be included in the potential solver. However, for 
the flow cases presented here, these terms had little effect on the 
solution, and, hence, for the results of this effort they were omitted 
from the governing equation. The next chapter presents and discusses 
details of the results obtained. Solutions are presented for an annular 
cascade consisting of eight bicircular arc blades with a maximum 
thickness of five percent. These blades are run for two geometric 
operating conditions: a high advance ratio case where the advancing 

speed is large compared to the rotational speed of the propeller tip, 
and a low advance ratio case where these two speeds are equal . For the 
high advance ratio case, the helical free-stream Mach number relative 
to the blade tip is M R = 0.8. For the low advance ratio case, results 
for two values of the free-stream Mach number are presented — a subsonic 
free-stream case of M R = 0.8, and a transonic case of M R = 1.1 which 
is representative of a turboprop. In general, the results are for zero 
angle of attack; however, a lifting case is presented for the low 
advance ratio and high free-stream condition. The nonlifting cases are 
compared with results from an Euler code (Ref. 6). A final chapter 
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delineates conclusions reached and recommends future efforts. Three 
appendixes are Included. The first two concern helical coordinate 
systems; the first of these presents the metric tensor for periodical 
coordinates used in the present work and the second presents 
nonperiodical coordinates used for flow about isolated blades. The 
remaining appendix contains a listing of the formulas for the derivative 
expressions relating stretched coordinates and uniformly spaced 
computational coordinates. 



II. PROPELLER CHARACTERISTICS AND NOMENCLATURE 

Some basic concepts will now be introduced related to common 

propeller characteristics and nomenclature. In many ways it is useful 

to think of the propeller as a highly twisted wing acted on by both 

lifting forces and drag forces. Lift is produced normal to the relative 

velocity at the surface. For small angles of incidence a., the lift 

is proportional to this angle. The angle of incidence equals the angle 

of attack a. minus the angle of zero lift a Q , as shown in 

Fig. 2.1. In a real fluid, drag is produced in a direction parallel to 

the relative velocity between the fluid and the body. The main 

distinction between a wing and propeller is that a propeller rotates and 

produces thrust while working its way through the fluid whereas a wing 

just advances while producing lift. The combined rotational and forward 

motion of the propeller results in highly twisted blades which are 

formed of radially varying profiles with the thickest sections nearest 

the hub. This increase in complexity admits additional parameters 

beyond those encountered in wing theory. 

The most important propeller parameter is the advance ratio which 

relates the forward velocity V to the rotational speed. The forward 

velocity is assumed to lie in the direction of the axis of rotation and 

will often be referred to as the axial velocity. The advance ratio J 

is commonly oefined as equaling the propeller's forward velocity divided 

by the product of its rate of revolution n times its diameter D: 

, V_ (2.1) 

J = nD • 
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FLIGHT DIRECTION 


FIGURE 2.1. - AIRFOIL SECTION SHOWING ANGLE OF ATTACK a. 
ANGLE OF INCIDENCE dj, AND ANGLE OF ZERO LIFT a Q . 
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Another definition of the advance ratio is frequently used. It is 
designated by X and equals the ratio formed by dividing V by the 
angular velocity of the propeller's tip 



where Q is the angular frequency of shaft rotation and R is the 
tip radius of the propeller. Convenience or mode of practice determines 
which definition is used since the two are related by 

J = irX. (2.3) 

In the present work, X will be used throughout. 

The re.ationship between the forward and rotational velocities for 
any propeller section is expressed in terms of the advance angle 13 
diagrammed and related to a helical curve in Fig. 2.2. The following 
relation connects the advance angle with the radius: 

tan 6 = ^ . (2.4) 

The slope of any helix is thus proportional to 1/r. The distance 
advanced during one revolution is termed the pitch, obviously the same 
for all radial elements of a rigid body. The surface swept out by the 
leading edge of a rotating and advancing propeller is a helical surface 
or helicoid and is called the advance helicoid. The angle of attack of 
any blade element is generally referenced to the advance angle 6; i.e., 
the angle of attack is measured from the advance helicoid. This is 
because a very thin blade constructed to lie in the advance helicoid 
produces no lift, and thus it gives a definition consistent with that of 
the planar wing. Blades are constructed such that the angle of attack 
varies radially. This variation in a can be interpreted as a radial 
variation in blade pitch since each radial element of the blade has a 
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geometrical pitch which, In general, differs from the pitch of the 

advance helicoid; the mean surface of a lifting, symmetric blade does 

not lie In a helical surface. This spanwlse variation of blade pitch 

introduces lift and allows the blade to be optimized in regard to its 

ratio of lift to drag. In elementary design the pitch of each section 

is chosen to give the optimum angle of attack for that section's 

profile. Figure 2.3 illustrates a typical span station showing the 

chord 8., angle of attack a, and the maximum thickness t v along 

m& x 

the chord. 

The performance of a propeller is measured in terms of its 

efficiency. The efficiency n is defined by 

THRUST x FREESTREAM VELOCITY (? ^ 

n " SHAFT POWER 

It is convenient to introduce dimensionless coefficients for thrust and 
power that have the form 

C T = and C P = ’ <2 ' 6) 

pn 0 pn D 

respectively, where T is the thrust, P is the shaft power, and p 
is fluid density. Then the efficiency can be simply expressed in terms 
of these coefficients and the advance ratio J in the following form: 



When stating the efficiency for a complete propeller system, it is 
necessary to denote the thrust, thrust coefficient and the corresponding 
efficiency in a manner indicating that the effects of components such as 
the spinner and nacelle are included. This is customarily done by using 
the subscript "net," such as in the following: 
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FIGURE 2.3. - AIRFOIL SECTION SHOWING CHORD L 
ANGLE OF ATTACK a, AND MAXIMUM BLADE THICK- 

NESS ^x- 
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n net 
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C T,net 
C P,net ’ 


III. HISTORICAL REVIEW 


3.1 Analytical Methods 

Investigations of flow about propellers have their roots in the 
early work of Rankine and Froude which was based on simple momentum 
theory for an incompressible flow. Simple application of momentum 
theory considers the momentum and energy of the fluid and treats the 
propeller as an actuator disk that imparts momentum to the fluid which 
passes through. This causes the formation of a slipstream of increased 
axial velocity. In performing this analysis (Refs. 7 to 8), it is 
useful to assume that the propeller system is composed of a large number 
of equally spaced blades producing a circumferentially uniform flow in 
which the velocity is constant across the disk but greater than that of 
the surrounding free-stream. A pressure change occurs across the disk 
which accounts for the momentum gain. By assuming irrotational flow and 
by applying Bernoulli's equation to streamlines on each side of the 
disk, this change in pressure can be related to the velocity rise in the 
final slipstream. The values of the pressure and velocity at pertinent 
locations are shown in Fig. 3.1 which illustrates this simple model. 

The value of the velocity increase, v, at the disk is found to be 
exactly half the increase, v 1 , existing in the final slipstream. 

Further results concerning propeller efficiency may be presented 
after considering the increase in axial momentum and rate of energy 
input. The use of momentum theory leads to limits on the maximum 
efficiency. Explicit relations are obtained between efficiency and both 
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ACTUATOR DISK 



FIGURE 3.1. - ACTUATOR DISK WITH UNIFORM CIRCUMFEREN 
TIAL FLOW. 
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power and thrust, which are customarily represented in terms of their 
coefficients, as a function of the most important propeller parameter, 
the advance ratio J. A more exact account of momentum theory, 
including the effect of slipstream rotation is given by Betz (Ref. 9). 
However, it should be mentioned that momentum theory fails to yield the 
detailed forces (loading) acting on each blade element. 

The attempt to determine the forces acting on the blades leads to 
the creation of the blade element theory. This theory approximates the 
aerodynamic forces acting on a blade by treating each radial element 
between r and r + dr as a two-dimensional wing section. The use of 
wing section theory provides the elemental lift and drag forces, which 
can be resolved in component form in the directions of thrust and 
torque, at each radial section. The values of the total thrust and 
torque can be determined by integrating the contributions of the radial 
elements. Corrections to account for the influence of neighboring 
elements can be added in an empirical manner. 

The application of blade element theory has been used in 
conjunction with momentum theory to provide improved results over 
solitary use of either method. This combined technique provides perhaps 
the best method for evaluating propellers with a large number of blades 
without actually calculating the entire three-dimensional velocity and 
pressure fields comprising the flow. However, for the normal situation 
involving propellers with few blades, the assumption of axial symmetric 
flow arising from the momentum viewpoint is weak since most of the flow 
does not pass near a blade. 

The evolution of propeller theory advanced significantly from the 
work of Goldstein (Ref. 10) which also stemmed from the successful 
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development o F the theory of the airfoil. This theory acknowledged the 
role of circulation in the generation of lift, a concept set forth by 
Kutta and Joukowski , and it also contained the formulation of Prandtl's 
model of flow past a finite wing of large aspect ratio. Prandtl's model 
introduced the concept of a continuously trailing vortex sheet which 
replaced the simpler model of two trailing tip vortex lines offered by 
Lanchester. c urthermore, this model allowed the actual airfoil to be 
replaced at its mean location by a bound vortex sheet across which a 
pressure difference can exist. For wings of large enough aspect ratio, 
the vortex sheet representing the wing can be approximated by a single 
bound vortex line lying in the plane of the wing normal to the 
approaching flow and with spanwise varying strength. Together the 
streamwise aligned vortices of the trailing vortex sheet and the lifting 
line vortex, representing the wing, induce velocities which superpose 
with the undisturbed stream velocity. The application of these 
propositions to the case of the propeller, including the modifications 
needed for three dimensional flows containing discontinuities, results 
in a complicated set of integral equations. 

Using the groundwork provided by early airfoil theory, Goldstein 

constructed a model which assumed that the trailing vortices emanating 

from the rear edge of the propeller forms a rigid vortex sheet of 

constant velocity, which closely approximates the structure of a helical 

sheet. The rigid vortex sheet was assumed to be embedded in an inviscid 

fluid with zero circulation. The flow far downstream of the propeller 

2 

was taken as satisfying Laplace's equation: V 4> = 0, where 

represents the velocity potential whose gradient is the fluid velocity. 
By aptly defining a helical coordinate along the screw surface in terms 
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of cylindrical polar coordinates, Goldstein was able to simplify the 
treatment of the boundary conditions. Satisfying single-valuedness and 
meeting continuity requirements, etc., led to a solution for the 
velocity potential far downstream of the propeller in terms of Bessel 
functions. Formulas were then derived from which the optimum load 
distribution along the blades could be calculated for any number of 
blades . 

Proceeding along the lines set by Goldstein, Reissner (Ref. 11 

to 12) derived and presented in more detail (Ref. 13) integral relations 

connecting circulation, thrust and torque along the span with the 

geometric and flow parameters such as blade angle distribution, number 

of blades and inflow angle. It is interesting to note that while the 

complete set ov equations presented was derived on the foundations 

developed from wing theory for incompressible flow, nevertheless, ideas 

introduced by H. Lamb (Ref. 14), for treating the magnetic field arising 

from a helical electric streamline, were drawn upon. Lamb presented the 

three-dimensional Biot-Savart integral in a manner which took advantage 

of helical symmetry and managed to reach solutions for the electric 

potential problem which were comprised of Bessel functions of the first 

kind for the inner radial field and of the second kind for the outer 

field. Reissner set out from this point, aided by the work of Kawada 

(Ref. 15), which was concerned with the Induced velocities resulting 

from a circulation of constant strength along a blade and which also 

made use of helical symmetry. Reissner based his mathematical 

2 

development on the Poisson differential equation, V 4> = q, for a 
velocity potential produced by an arbitrary source distribution. He 
proceeded toward a more general integral equation, which differed from 
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the Biot-Savart integral and the series solution obtained by Goldstein, 
by replacing the blades with lifting lines with infinitely thin trailing 
vortex sheets. The vortex sheets were represented by the proper 
distribution of sources and sinks entering on the right-hand side of 
Poisson's equation. Helical symmetry was invoked to reduce the number 
of independent coordinates to two, with the new independent coordinates 
being a radial coordinate and a helical coordinate defined in the manner 
of Goldstein and Kawada from cylindrical coordinates. The ability to 
arrive at this two-dimensional equation for the velocity potential 
relied on the vortex sheets being arranged so as to lie on multiple 
sheets of helical symmetry with Infinite extent downstream. Further 
refinements, which removed this symmetry condition but improved the 
theory, were pointed out by Reissner. These included taking into 
account the radial contraction of the flow on passing through the 
propeller system (slip-stream contraction) and the displacement of the 
helical vortex sheets due to their mutual interaction. These 
Improvements were considered, however, to be of a secondary nature. 

Shortly before the end of World War II, Theodorsen published a 
series of technical papers (Refs. 16 to 19) dealing with both single and 
dual rotation (counter-rotating) propellers. These reports formed the 
basis of his subsequent book Theory of Propellers (Ref. 20), which 
presented experimental measurements as well as the theoretical 
treatment. 

Theodorsen extended the lightly-loaded theory presented by 
Goldstein in his 1929 paper to that of heavy loading by redefining 
certain key propeller parameters. Goldstein had acknowledged that a 
better approximation would be gained, valid for moderate loads, by doing 
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such. This modification took Into account the Increase In advance ratio 
and slipstream contraction that can occur to a wake surface as It 
travels downstream; see Fig. 3.2. Expressions were derived for the 
thrust, torque, and efficiency in terms of these parameters. Included 
with these parameters were the circulation coefficient connecting the 
circulation at each blade element to the advance ''•atio and the mass 
coefficient parameter which represented an average of the circulation 
coefficient over the cross-section normal to the flow. It is important 
to note that these parameters referred to the ultimate wake and not to 
the region immediately behind the blades. Theodorsen reasoned that the 
distribution of circulation existing on the propeller blades is, to the 
first order, the same as the resulting distribution of potential 
difference existing in the wake if the radius of the wake is stretched 
at each axial location to equal the radius of the propeller blades. The 
optimum circulation distribution for heavy loading is thus identical to 
that of light loading provided the helix angles in the ultimate wake are 
the same. Thus, the key relations for propeller theory, as developed 
along the lines of Goldstein, were extended to heavy loading and 
expressed compactly in terms of the important mass coefficient 
parameter. 

In addition to the theoretical development, Theodorsen carried out 
extensive detailed measurements. Here he relied on the analogy between 
electrical theory and fluid mechanics for potential flow to conduct 
electrolytic-tank measurements in place of conventional wind-tunnel 
tests. By placing a helical insulating surface into an electrolytic 
tank with a uniformly imposed electric field along it axis, an electric 
field entirely descriptive of the discontinuous flow behind a propeller 
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FIGURE 3.2. - WAKE PROPAGATING DOWNSTREAM CONTRACTS 
CORRESPONDING TO AN INCREASE IN PITCH. 
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was produced. The effect of advance ratio and number of blades on the 
mass coefficient was readily determined by varying the helix shape and 
geometry of the tank. This provided the experimental results needed for 
comparison with theoretical results and extended the results beyond the 
known range to include higher advance ratios. The general procedure is 
covered in some detail for applying the theory to the general design 
problem. The interested reader is referred to Theodorsen's book. 

Still within the structure of linearized theory, Busemann 
(Ref. 21), and Davidson (Ref. 22) attempted to include compressibility 
effects by generalizing the incompressible propeller theory to 
supersonic tip speeds. These two investigations progressed from origins 
similar to those of Goldstein and Reissner. Busemann advanced physical 
arguments to explain or predict changes in the flow for higher Mach 
numbers; Davidson provided detailed calculations as well. The most 
notable characteristic of flows at Mach numbers approaching unity is the 
mixed nature of the flow fields. For mixed flows the elliptic character 
of the solutions for subsonic regions, which are essentially 
incompressible, are joined by the hyperbolic character representative of 
the supersonic regions where compressibility effects cannot be ignored. 
The existence of supersonic as well as subsonic regions of the flow was 
recognized and tied to the solutions which again took the form of 
Fourier-Bessel functions. In these cases the solutions obtained 
satisfied the steady linearized small disturbance equations. Busemann's 
equations were restricted to two dimensions involving polar coordinates. 
Davidson developed a more general disturbance equation which was 
expressed in both cylindrical and helical coordinates. Although the 
equations contained no nonlinear terms, they exhibited the mixed nature 
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of the flow field which progressed with increasing radius from an 
entirely subsonic (subcritical ) flow to a flow embedded with pockets of 
supersonic (supercritical) flow to, finally, totally supersonic regions 
at the tip. In contrast, the fixed wing problem results in essentially 
only one of the above cases depending on the approaching Mach number. 

It is largely a result of the radial variation of the helical velocity 
that makes tha propeller problem more complex than that of the wing. 

This is evidenced by the presence of the spanwi se-type mixed flow 
feature even with linear theory. 

Even within the limitations of linearization, the effects of 
numerous geometric parameters must be taken into account for a more 
rigorous treatment of propeller flow. The effects of thickness, camber, 
angle of attack, sweep, offset, blade interference, and tip relief must 
be dealt witn to a sufficient degree. A number of more recent 
investigations has addressed different aspects of these problems. In 
particular, an analytical effort by Hanson (Ref. 23) to include these 
items resulted in a fairly general description of the flow field. He 
constructed an analysis for prop-fans based on the acoustic methods of 
Goldstein (Ref. 24) which relate the density disturbance at any point 
in a medium enclosed by impermeable walls to the arbitrary motion of 
unsteady sources. This is accomplished via an integral formulation 
employing Green's functions and provides an integral equation for the 
near and far field density disturbance and one for the unsteady velocity 
potential. The result for the disturbance pressure, or acceleration 
potential as it is often called, is valid for any Mach number and 
planform but presumes that the surface pressures are known. The 
formulation for the derivation of the velocity potential contains only 



31 


linear source terms and is valid for any steady or unsteady helically 
convected source. The velocity potential itself must be obtained by 
inversion methods used for wing analysis. Two key elements were needed 
to achieve these results. First, the exact boundary conditions at the 
blade surfaces were actually applied at mean reference surfaces 
approximating the blade positions, thus greatly simpMfying the 
mathematics. Secondly, it was recognized that a nelicoidal surface 
theory could be used as a natural extension of wing theory with the 
observation that each blade element rotating in forward flight sweeps 
out a helix. Thus, a coordinate system was established that was formed 
by a dimensional coordinate at constant radius in the advance helicoid, 
a radial coordinate, and another coordinate at constant radius lying in 
the helicoid normal to the advance helicoid. Blade coordinates are then 
transferred to the advancing helicoid at locations neighboring the upper 
and lower blade surfaces as in thin airfoil theory. In this manner the 
ideas of wing theory were adapted to propeller analysis. This analysis 
then gives aerodynamic and acoustic information pertaining to steady 
performance, unstalled flutter, and noise radiation. 

This discussion concludes the section on analysis methods; it is 
necessarily incomplete, but presents a review introducing many of the 
important concepts developed earlier in propeller theory. It should be 
stressed that the above methods neglected many physical phenomena in 
order to simplify the mathematical models and provide solutions. 

Indeed, their success hinged mainly on construction of simplified 
models. These models retain most of the significant physics, and which 
could be solved and lead to improved designs. It is testimony to the 
insight and rasourceful ness of these early researchers that such large 


strides could be made. 
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3.2 Nonlinear Flows with Shocks 

This section discusses the development of CFD methods applied to 
high-speed flows where compressibility effects are important and in which 
shocks may be present. Generally the effects of viscosity and heat 
conduction will be considered negligible, though some comments will be 
made concerning the former. The development of successful CFD methods 
has relied on an understanding of the physical characteristics and the 
mathematical properties of the partial differential equations associated 
with these flows. The physical significance and the mathematical 
behavior of the governing equations for inviscid flows of this type will 
be reviewed in order to establish the consequences of transitioning 
from partial differential equations to finite difference equations. 

Difference equations can be thought of as approximating the partial 
differential equations (PDE) on a general lattice of mesh points. In 
this manner, the continuous functions present in the PDE will be 
replaced with finite difference functions that are discrete and defined 
only on the lattice and only at discrete intervals of time. The finite 
difference equations will be considered consistent with the correspond- 
ing PDE if it can be demonstrated that any difference between the PDE 
and the finite difference representation vanishes as the general 
lattice of mesh points is refined. This can be accomplished if a 
substitution of an assumed solution to the PDE in the form of a Taylor 
series expansion into the finite difference equation produces agreement 
to within a set of higher order terms which tend to zero as the mesh 
spacing and time interval simultaneously are allowed to approach zero. 

The higher order terms are called the truncation error terms and 
represent one measure of solution accuracy. The accuracy cannot be 
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predicted beforehand since in each truncation error term, a derivative 
factor of differential order corresponding to the order of the term is 
also present, but unknown. For cases where large gradients in the 
solution exist, these derivatives can be large enough to offset the 
smallness of the mesh. Estimates of accuracy can be made once a 
solution is achieved by examining the truncation terms or by repeating 
the calculation for several levels of mesh density. The consistency 
condition is an important requirement for a difference scheme. If a 
difference scheme converges, then it will converge to the solution of 
the corresponding differential problem provided that they are 
consistent with each other. 

In addition to consistency, one would prefer that the solutions for 
a set of finite difference equations possess certain properties desired 
similarly, though not necessarily guaranteed, for the PDE. One hopes 
that a solution exists that is unique and stable. In this regard a 
problem can be considered to be stated properly if it has exactly one 
stable solution satisfying both the governing equations and the 
auxiliary conditions. The auxiliary conditions are usually boundary or 
initial conditions but may include both or other types of constraints. 
The existence and uniqueness requirements are dependent on the stated 
conditions being consistent and sufficient to determine a unique 
solution. The stability property manifests itself in different 
contexts. A stable solution of a mathematical physics problem implies 
that, for small changes in a given condition, there will be 
correspondingly small changes in the solution. The question of 
stability is very important since inevitably the process of carrying out 
a solution introduces some error through truncation or the inexact 
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specification of some condition. It is essential that the effect of 
these errors remain small, causing only small inaccuracies in the 
solution. 

The solution of a problem using a finite difference grid carries 
with it inherent inaccuracies. Since information is utilized only at 
discrete points, small wavelength components are not resolved 
sufficiently, and those shorter than the length of a grid interval not 
at all. The question of stability in regard to finite difference 
equations can be addressed in terms of the growth in amplitude of 
frequency harmonics of the solution, and it is subsequently related to 
the ratio of the coarseness of the physical grid and the time step 
increment, in say an initial value problem. The iterative updating of 
the solution from a given time level to a new level At later is called 
the elemental marching step. Depending on whether the numerical scheme 
is of the implicit or explicit type or a mixture of tnese two types, 
certain frequency harmonics can experience unlimited growth in amplitude 
as the elemental marching steps are performed. The numerical scheme is 
said to be unstable if any frequency dominates the solution in such a 
physically unwarranted fashion. Restricting the time step to conform to 
the Courant-Friedrichs-Lewy (CFL) (Ref. 25) condition generally assures 
that such components will not grow exponentially, as long as additional 
destabilizing influences do not exist beyond those associated with the 
elemental marching step. Certain methods of enforcing the boundary 
conditions may lead to instability regardless of satisfaction of the CFL 
condition. Exact conditions for stability are not attainable for most 
time dependent and nonlinear problems. Recourse is often made to linear 
stability criteria obtained by neglecting nonlinear effects. This is 
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helpful in determining guidelines on stability bounds in many instances. 
The existence of stability and satisfaction of the consistency condition 
are sufficient to insure that the finite difference solution converges 
to the solution of a linear differential equation as the difference 
intervals tend to zero (Lax's Equivalence Theorem). Unfortunately, a 
similar statement cannot be made for nonlinear problems. The properties 
of existence, uniqueness, and stability are very important, but the 
problem of determining whether a problem is stated correctly is also 
very difficult. Some of these questions will be addressed further. 

The nonlinearity of the inviscid transonic flow equations is the 
main hindrance In the development of analytical methods capable of 
solving these partial differential equations. Few instances arise where 
this difficulty can be avoided. Use of the hodograph plane transforms 
the basic equations into a linear form by transposing the roles of the 
independent variables. However, the hodograph transformation is 
restricted to two-dimensional flows, and even then, a trade-off exists 
in that the boundary conditions are generally more complicated though 
the equations are linearized. Thus, in most cases, the transonic flow 
equations wil 1 be encountered in a nonlinear form. In addition, 
shock-free inviscid flow solutions are rare, and if they exist are 
likely to appear as isolated solutions which are close to solutions 
containing shocks. This section will discuss some of these aspects and 
how they enter Into consideration when devising finite difference 
methods applicable to transonic flows. With respect to numerical 
efforts, only a few cases will be considered which have a bearing on 
the present wo^k. 
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Much of the progress achieved in the development of numerical 
methods for steady flows has been found to extend directly to unsteady 
applications. In this regard, steady flows are often solved as 
converged solutions of unsteady problems, asymptotically reaching 
steady state through real or pseudo time. Mention will be made of both 
unsteady as well as steady inviscid flows. Since this work solves a 
hyperbolic type equation due to the introduction of time as another 
independent variable, the unsteady perspective will be stressed. The 
efficient and accurate solution of a given flow, whether steady or 
unsteady, requires proper numerical simulation of the dominant physical 
phenomena. Hence, with the understanding that the proper numerical 
approach is more fully appreciated if the physical aspects of the 
problem are first understood, additional features pertaining to 
transonic flows should be stated beyond merely recognizing their 
nonlinearity. Transonic flows are characterized by the following 
traits: sensitivity to small perturbations; weak decay of these 

fluctuations as they propagate laterally; slow global adjustment to 
local disturbances; mixed nature with supersonic regions embedded within 
subsonic regions; and likely presence of surfaces of discontinuity such 
as shocks. Simple and accurate treatment of shocks is of primary 
importance in providing efficient numerical schemes for transonic flow. 
Therefore, a brief description of their nature will be given along with 
a quick review of how they are incorporated into the numerical modeling. 

A simplified picture of shock waves can be provided if viscous and 
heat conduction effects are neglected while considering the response of 
a flow to a disturbance. The placement of a thin body in a uniform 
stream which approaches it at high subsonic or supei sonic speeds 
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produces a discontinuity in the pressure field. This discontinuity 
results because high pressure pulses propagate more rapidly than low 
pressure ones and shock waves invariably form by coalescing pressure 
waves originating from any finite disturbance. The pressure waves must 
be compression waves in order for a shock to form, never expansion 
waves. In actuality these discontinuities are smeared out by physical 
dissipative effects which are absent when viscous and heat conduction 
effects are discounted; only in the limit of zero viscosity do shocks 
appear as discontinuities. However, the mathematical description of a 
shock is consiaerably simplified by neglecting these dissipative 
mechanisms and the approximation to reality is quite good in most cases. 

It can be shown for unsteady plane flow that any disturbance that 
compresses the medium will eventually lead to a breakdown in the gas, 
i.e., a discontinuity will develop when a simple plane wave produced by 
a compression force propagates into a gas at rest. An example 
consisting of a limiting form of Burger's equation will illustrate the 
formation of a discontinuity from an initially smooth wave. It will 
also serve the purpose of introducing some concepts related to nonlinear 
partial differential equations that are basic to their understanding. 

In fact, many of the fundamental investigations concerning numerical 
computations deal with this particular equation, and thus it serves as 
a natural springboard to pass to more complicated equations by providing 
valuable core material underlying this area. 

Burger's equation is a widely studied nonlinear equation modeling 
both convective and diffusive processes and is of the form 
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where u, x, t, and u can be taken as the nondimensional velocity, 
distance, time, and diffusion coefficient, respectively. Actually 
Eq. (3.1) is more accurately classified as a quasi-linear equation since 
It Is linear with respect to the highest order derivative terms. We 
wish to examine the special case of u = 0; no diffusion occurs and 


Eq. (3.1) reduces to the homogeneous form 
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ax 
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(3.2) 


Specifying an initial condition 

u( x , t ) t=0 = u(x ,0) = sin(x) (3.3) 

defines an initial value problem with a smooth initial curve. The 
solution to this problem is given in implicit form by 

u - s i n ( x - ut) = 0 (3.4) 

which is readily verified by its substitution into Eqs. (3.2) and (3.3). 
This solution is unique but exists only within a limited range of the 
Initial curve. The solution is sketched in Fig. 3.3 for several time 
intervals. The waveform of the solution can be seen to progressively 
distort from the smooth shape at time t = 0 to a discontinuous shape 
at time t = 2 in the manner described above. This demonstrates an 
Important property of the solution to nonlinear equations: they can 

develop discontinuities even with smooth initial data; this is in 
contrast to the linear case where discontinuities can arise only if the 
initial data are discontinuous. This indicates the necessity for 
admitting possible discontinuous solutions when solving nonlinear 
problems. More will be said later concerning discontinuous solutions 
to the equations of gas dynamics. 


RE 3.3 
OOTH Wi 
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We return to the examination of Burger's equation by connecting it 
to the theory of partial differential equations. This will provide 
access to some valuable geometric and analytic insight. To carry this 
out, observe that equation (3.2) is a quasi-linear differential equation 
in two independent variables, x and y, of a general form 

P fx + Q fy = R (3 ' 5) 

where P, Q, and R may be functions of x, y, and u only; 
furthermore, all variables and their first derivatives are continuous 
and P 2 + Q 2 * 0. In the special case presented above, P = u, Q = 1, 

and R = 0, and where in Eq. (3.2), t has taken on the role of the 

variable y that appears in the more general Eq. (3.5). 

A general theory exists for equations of this type which can be 
drawn on in order to provide analytical results (Refs. 26 to 27). The 
solution to Eq. (3.5) can be interpreted in terms of an integral surface 
defined through the function 

F(x,y,u) = C (3.6) 

where C is a constant and u satisfies Eq. (3.5); such a surface can 
be visualized geometrically in an x,y,u-space. The most general 
solution to Eq. (3.5) is a function G(u^,u 2 ) where and u 2 each 

satisfy the equation and correspond to separate values of the constant 
C, say C 1 and C 2> appearing in Eq. (3.6). Thus u 1 and u 2 are 

two distinct integral surfaces. The function F can be introduced into 
Eq. (3.5) resulting in an equation in which u can be considered 
independent along with x and y. The partial derivatives of F with 
respect to both x and y separately are each equal to zero as can be 
seen immediately from Eq. (3.6). These can each be used to give 
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expressions for one of the partial derivatives appearing in Eq. (3.5) 
which then assumes the following symmetrical form 

D 3F n 3F D 3F n 

P ax +Q 3y tR 3u’° <3 ' 7> 

provided that 

This can be thought of as an expansion of the scalar product formed from 
vectors having components P, 0, R and 3F/9x, 3F/3y, 3F/3u in some 
appropriate orthogonal Cartesian system. Geometrically this has 
significance in that it shows that the vector with components P, 0, R 
is perpendicular to the normal to the integral surface at any point. 
Thus, the vector lies in the plane tangent to the surface at any point 
and its direction specifies a fundamentally important direction on the 
surface called the characteristic direction. This direction is unique 
at any point but, in general, differs with location. Starting at any 
point in the surface and following the characteristic direction along 
the surface results in a curve being traced out called the 
characteristic curve; hence, a characteristic curve is everywhere 
tangent to the characteristic direction. The projection of the 
characteristic curves onto the x,y-plane defines the physical 
characteristics. 

An equivalent way of expressing Eq. (3.5) can be given that 
promotes its visualization geometrically by introducing certain 
parameters. This will aid in the construction of solutions which pass 
through a given curve. Extension of these concepts to higher dimensions 
provides an approach to constructing solutions for higher order problems 
by the method of iteration. From the theory of quasi-linear 
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differential equations, the characteristic curves are defined in terms 
of the following set of first order ordinary differential equations 

§§ = PCx.y.u) 

jj* = Q<x,y,u> (3.8) 

^ = R(x,y,u) 

where s is a continuous parameter varying monotoni cal ly along these 
curves. A two-parameter family of integral curves satisfies Eq. (3.8) 
since there are two independent differential equations in this set, 
i.e., the set cf equations depends on s and two additional constants. 
It is convenient to introduce new variables such that one is s which 
varies along each characteristic and the second variable x is 
introduced as a parameter which is constant along any characteristic 
curve but varies from curve to curve. Together, the curves form a one- 
parameter family of solutions sweeping out the solution surface defined 
in Eq. (3.6) for a given value of C. Since the solution of the 
ordinary differential system admits a two-parameter family of solutions, 
an arbitrary function must be introduced to define a particular one- 
parameter family subset which satisfies Eq. (3.5). This arbitrary 
function relates the two arbitrary constants that arise from the 
sol ut ion of Eq . (3.8). 

Specification of this function gives rise to an initial value 
problem for which the requirement that the solution surface pass through 
a particular curve establishes an initial value problem and fixes one 
parameter. The original independent variables can be defined in terms 
of the new variables as 


43 


X = x( s ,-c) 
y = y(s,x) 
u = u(s , t) 

such that x, y, and u satisfy the initial condition that 

x Q = x(0,x) 
y Q = y<0,x> 

= u(0 ,x) . 
o 


(3.9) 


(3.10) 


Here the value of s along the initial curve <€ is arbitrarily chosen 
as zero. It should be mentioned that the data curve x , y , u needs 
only to be continuously differentiable. 

Solutions to the above system exist in some neighborhood of the 
given initial curve if the above transformation has the property that 
the Jacobian 

J = xv - y e x (3.11) 

does not equal zero. This guarantees that the transformation is 
invertible such that s and x can be expressed in terms of x and 
y. In this case the solution u(s,t) of Eq. (3.9) can be expressed as 
u(x,y) and is generated by the characteristic curves which must lie 
entirely in the surface. Should the Jacobian vanish along < <f, then 
either # is a characteristic curve and an infinite set of solutions 
exist or it is a noncharacteristic curve and no continuously 
differential solution exists. This statement applies to any other curve 
displaced from the initial curve. In general, should any two solutions 
u 1 and u~ pass through the same curve, then it must be a 
characteristic curve. The significance of these remarks will become 
clear when requirements are discussed concerning the proper treatment of 
shocks as discontinuities. 
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By returning to the example of Burger's equation, the parametric 
approach can be applied directly giving information on the solution 
domain. Assume that a solution surface exists satisfying the 
homogeneous Burger's Eq. (3.2) which is generated by a one-parameter 
family of curves passing through the initial curve defined by Eq. (3.3). 
The first two expressions in Eq. (3.8) define the family of 
characteristics and the last one defines the solution surface in terms 
of the parameter which is constant along each curve. This can be can 
seen by carrying out the integration of Eq. (3.8) by recalling that 
P=u, Q * 1, y = t and, since it is homogeneous, R = 0. The solution 
which satisfies the equation and the condition that it initially be a 
sine curve is given parametrically as 

x = s sin x + x 

t = s (3.12) 

u = sin x 

(Ref. 28). The Jacobian of the transformation is 

X s y x " y s X x = _(1 + s cos *>• (3.13) 

It can be seen that for certain values of s and x the Jacobian 

will be zero. The locations where it vanishes are of interest, and the 
projection of these locations on the physical plane is termed the 
envelope. Such an envelope denotes positions where projections of the 
characteristics cross which, in turn, implies that the solution is no 
longer single-valued. The location of the envelope Tor this example can 
be determined by introducing the variables x and t into Eq. (3.13) 
and then setting this equation for the Jacobian equal to zero. The 
resulting curve for the envelope is then given by the relationship: 


x 


(3.14) 
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The solution for the envelope Is sketched In Fig. 3.4. Since the value 
of R Is zero for this problem, the third expression in Eq. (3.8) shows 
that the solution is constant along a characteristic curve which is 
verified by the last relation in Eq. (3.12). In the region to the right 
of the curves AB and AC in Fig. 3.4, characteristic curves cross and 
Eq. (3.14) is multi-valued; this corresponds to the dotted portion of 
the curve at t = 2 shown in Fig. 3.3. The discontinuity arising in 
the solution is a result of the crossing of characteristic curves, each 
associated with a constant value of u. Since the value of u is 
different for characteristics intersecting at the discontinuity, there 
is a resulting jump in the value of u associated with these 
conflicting characteristics. This exhibits a general feature of 
nonlinear equations. A discontinuity, which may be present in a 
nonlinear flow, does not lie along a characteristic curve. This is in 
contrast to a linear problem where discontinuities are present only if 
the initial data are discontinuous, in which case they are propagated 
along the characteristics. More specifically, it can be stated that 
shocks do not lie on characteristic curves but rather separate the flow 
Into distinct regions differing by jumps in certain physical quantities. 
The jump In the physical quantities must satisfy certain physical 
criteria. It is exactly these jump criteria which serve to sift out 
the physically correct solution from the other mathematical solutions. 

Because of the presence of these discontinuities, it is necessary 
to relax the stipulation that a solution be continuously differentiable. 
This fact stems from consideration of integral conservation laws which 
equate the rate of change in a quantity in a region to its flux at the 
boundary of the region. We can expand the definition of u to that of 
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a generalized density of a quantity and denote its flux as f. The 
integral form cf the conservation law for the quantity is 


d_ 

dt 


/ u dx = - I f • 

Jq •'s 


ndS 


(3.15) 


for a general domain of arbitrary dimension with bounding surface S 
possessing an outward normal n and where dx signifies the 
differential ''volume" element. Here, we follow the presentation of Lax 
(Ref. 29) who omits the vector symbols in the scalar product. The 
conservation equation states that the rate of increase of the quantity 
in the region equals the negative of its outflow through the bounding 
surface. This holds for quantities like mass, momentum, and energy but 
not, in the general case, for entropy or internal energy for example. 

Equation (3.15) is a general relationship to which the divergence 
theorem can be applied to the right-hand side, for sufficiently well 
behaved functions, which transforms it to 


<u t + div f) dx = 0 (3.16) 

with the time derivative moved inside the integral and denoted as a 
subscript. For u and f which possess continuous partial 
derivatives, the derivative form of the conservation law follows 



ut + di v f = 0. (3.17) 

This is the so-called divergence form of the differential conservation 
laws. The requirement that the derivatives be continuously 
differentiable implies that discontinuities such as shocks must be 
treated as internal boundaries separating the total domain into regions 
which by themselves possess smooth solutions; Eq. (3.17) then can be 
satisfied by these piecewise continuous solutions. Lax calls the 
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solution for the total region, including the discontinuity which 
satisfies the integral law, a generalized solution so as to distinguish 
it from a regular solution which must be differentiable. 

Recall the example of the sine wave given above. The solution near 
the initial curve begins as a generalized solution, which is regular, 
but evolves into a generalized solution that is no longer continuous. 

The concept of a discontinuous solution is fully justified provided that 
it is physically meaningful. However, discontinuous solutions which 
satisfy the differential conservation laws may not be unique. The 
criterion that is used to distinguish the physically relevant solution 
is that the jump conditions existing across the discontinuity must 
satisfy the physics of the problem. 

When shocKS are treated as discontinuities, the differential 
equations must be replaced with equivalent upstream and downstream 
boundary (or "jump") conditions. These are known as the Rankine- 
Hugoniot relations. They can be obtained from the conservation laws, 
the Euler equations, by applying these laws separately to a control 
volume which contains the shock and which moves at the shock velocity. 
These relations specify the changes in the thermodynamic properties 
such as temperature, pressure, and density occurring across a shock 
without considering the detailed physics within the shock. In this way 
physically meaningful solutions are obtained. The Ranki ne-Hugoniot 
relations, which conserve mass, momentum, and total energy, were 
actually preceded by a set of analogous equations in which entropy was 
conserved rather than energy. Since entropy is not conserved across a 
shock, these analogous relations are normally discarded in favor of the 
Ranki ne-Hugoniot relations. Since shocks are the only means of entropy 
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production in inviscid flows, either set of conservation equations would 
be valid in a flow without shocks. 

Some comments should be made concerning the potential equation 
derived from Euler equations by requiring the flow to be irrotational . 

By Crocco's theorem for steady flows, an irrotational flow is isentropic 
if the flow has constant total enthalpy; however, the flow may be 
rotational and still be isentropic. In this sense, isentropic flow is 
consistent with irrotational flow in that it imposes no additional 
restrictions on steady flows with uniform stagnation enthalpy. Thus 
potential calculations generally assume an isentropic jump, and 
consequently either conservation of mass, momentum, or energy must be 
sacrificed to prevent an overly determined system. The choice of 
conserving mass and energy is usually made, based mainly on the fact 
that they are both scalar quantities and that, for steady flows, 
integration of the energy equation implies constant enthalpy along 
streamlines. Klopher and Nixon (Ref. 30) review the three possible 
pairings for the conservation quantities and discuss the likelihood of 
conserving mass, momentum, and energy by allowing entropy changes to 
occur across the shock. 

Conservation laws can also be expressed for a system of equations 
and take the differential form 

+ div fj = 0 j = 1 ,2,3,.. .,n (3.18) 

dt * 

for n unknown quantities and their respective fluxes f J , where each 
flux is, in general, a nonlinear function of all the u J . This can be 
written in an equivalent form where the divergence operator is replaced 
by differential operators appearing as an implied double summation over 
both repeated dummy indices i and k as such: 
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(3.19) 


The index i ranges up to the number of independent variables not 
including time; the index k ranges over the number of n unknowns. 
This notation is closely connected to matrix notation. This can be seen 
by identifying u as well as the partial derivative factors in the 
terms with matrices. Let W be a column vector of the unknowns. Then 
the following matrices 
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lead to the compact matrix representation 

W t + A i W x i = 0 (3.21) 

for a system of n equations with the implied sum over each of the i 
independent variables other than time; the matrix A. is an n by n 
matrix. This notation is quite prevalent in the literature. The reason 
for introducing the system of quasi-linear equations is that the Euler 
equations are such a set of equations and the equation that will be 
solved herein is a limiting case of them. 

The Euler equations can easily be written in vector form for two- 


dimensional unsteady flow with velocity components u and v in the 
respective directions of x and y as 
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W t + F x + G y = ° (3.22) 

where 


W = 

and p is the density. Here, e is the specific internal energy and h 

2 2 2 

is the specific enthalpy for an ideal gas and q = u + v . The 
pressure is given by the equation of state p = pRT for a perfect gas at 
temperature T, where R is the universal gas constant. Obviously F 
and G take on the role of fl for i = 1 and i = 2, corresponding 
to the flux in the x and y directions, respectively. The attached 
subscripts denote partial differentiation with respect to the indicated 
variable. For isentropic flows, the pressure is related solely to the 
density through the relationship pp - ^ = constant where & is the 
ratio of specific heats. The Euler equations then simplify such that 
W, F, and G reduce to only their first three elements. 

The jump conditions relating the upstream and downstream states 
adjacent to a shock are given by Lax (Ref. 29) as 

• n (3.23) 

for a smooth shock surface moving with speed v s> The notation is the 
same as that used for Eq. (3.15); however, the brackets indicate the 
jump values of the enclosed quantities occurring across the shock. The 
value of the shock speed v $ is limited above by the speed of sound 
preceding the shock, and bounded below by the speed of sound on the high 





52 


pressure side. This guarantees that the shock will experience an 
entropy increase rather than a decrease, a feature that the Rankine- 
Hugoniot relations fail to discriminate. Since the initial conditions 
are not sufficient in themselves to determine a unique solution, these 
are the supplemental conditions needed to correctly specify a shock. 
Failure to use either the correct conservation equations or the 
conservation form leads either to incorrect shock speeds or strengths. 

The process of calculating the position of a shock and the 
subsequent application of the Rankine-Hugoniot relations is known as 
shock fitting. The use of shock fitting provides an accurate method of 
handling shocks and it also provides sharp resolution since they are 
treated as discontinuities. Shock fitting procedures were used as early 
as the year 1948 by Emmons (Ref. 31) who used a somewhat different 
approach in establishing shock jump conditions for flow over airfoils. 
The direct application of these fitting techniques is handicapped since 
the positions describing a shock surface must be calculated, often 
resulting in lengthy trial and error procedures. This is particularly 
troublesome since, in general, a shock moves through the medium along 
an undetermined path. Numerical schemes which fit shocks frequently 
exploit the theory of characteristics since certain quantities are 
constant, or nearly so, along the characteristic directions. 

The notion of characteristics which was mentioned above in regard 
to first order equations plays an important role in the development and 
understanding of numerical methods for systems of hyperbolic equations. 
To see how the theory of characteristics enters, it will be simply 
stated that the matrix for the system of equations in Eq. (3.21) can be 
transformed by multiplying it by a linear matrix which results in an 
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equivalent set of equations which have an interesting property. This 
property is that each of the resulting equations contain differentiation 
in a single direction. This is called the normal form for a system of 
hyperbolic equations. The directions of differentiation are the 
familiar characteristic directions which were presented above in the 
more easily understood form of first order equations with two 
independent variables. Nevertheless, the characteristics generalize in 
the case of higher order equations or systems of equations with more 
than two independent variables; the main difference is the obvious 
increase in diversity and complexity exhibited when more variables are 
Involved. Higher order equations are often treated as systems of first 
order equations by introducing extra dependent variables. 

For the unsteady one-dimensional Euler equations, three 
characteristic directions emanate from each point in the fluid. Along 
these characteristic curves information propagates, either in the form 
of waves or fluid motion. Sound waves flow along forward and rearward 
characteristics while particle paths follow a characteristic upon which 
the entropy is constant and which lies between the other two. For 
isentropic flow, certain quantities known as Riemann invariants can be 
constructed from the fluid variables such that one is constant for the 
forward characteristic and another for the rearward direction. These 
directions are sketched in Fig. 3.5. Signals reaching a point N, which 
just happens to coincide with a lattice point in the flow, have paths 
labeled \ + and representing the forward and rearward 

characteristics, respectively. In general, the characteristics do not 
follow paths connecting grid points. For noni sentropi c flow, quantities 
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FIGURE 3.5. - DOMAIN OF DEPENDENCE FOR POINT N LIES 
BETWEEN CHARACTERISTICS. 
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analogous to the Riemann invariants can be defined which are 
approximately constant, at least for small time periods. 

The fact that initial data propagated along characteristic 
directions determine the nearby solution to such hyperbolic problems 
leads to the important concept of domain of dependence. Data outside 
the characteristic rays shown in Fig. 3.5 cannot influence the solution 
at point N. Only data within the region influences the solution. This 
region is called the domain of dependence. Similarly, the 
characteristics leaving a point at time t will define a domain 
influenced by that point later on. It is important for stability 
reasons that these domains be respected while carrying out the finite 
difference calculations. The time-step must be adjusted to the spatial 
grid interval so that the analytical domain of dependence is contained 
within the computational domain of dependence. This restricts the 
maximum time-step on the basis of qualitative physical arguments. For 
example, inspection of Fig. 3.6 reveals that information from the three 
circled grid points at time level t would suffice to provide all 
necessary data for point N at time t + At. The computational domain, 
which is represented as straight dashed lines, contains the physical 
domain of dependence framed by the solid lines. For double the time- 
step, say in the absence of information at time t, then five grid 
points of data at time t - At would contain all the information 
influencing point N. The region between the solid and dashed lines on 
each side of the domain of dependence represents information which is 
not needed at point N but which it sees nevertheless. This transfer of 
spurious data will occur as long as the characteristics do not connect 
the grid points; in general they do not and, consequently, loss of 
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accuracy results. The domain of dependence principle places a. 
qualitative restriction on the time-step. This is of special importance 
to explicit numerical schemes. Implicit schemes have the advantage that 
time-step limitations imposed for stability may be relieved or removed 
altogether. In any case, stability requirements are usually not 
available for a nonlinear problem but hopefully estimates can be 
acquired from a linear counterpart or numerical experiment. 

Further inspection of the domain of dependence is required in the 
instance of embedded shocks. To illustrate this case, Fig. 3.7 presents 
a sketch similar to those just reviewed but which now includes a shock 
path in the computational plane. Assuming the flow to be supersonic on 
the left (the upstream side) of the shock and subsonic on the right, two 
sets of characteristics can be identified with points N and N+l residing 
in the supersonic and subsonic zones, respectively. In the supersonic 
zone, it is clear that all soundwaves travel downstream, and hence, both 
X-characteri sties have positive slopes; i.e., they both propagate 
downstream. For the subsonic zone, one characteristic points upstream 
and one downstream. Information affecting a supersonic point must come 
from the upstream direction, commonly referred to as the upwind 
direction. Subsonic points are influenced by signals from each 
direction. Finite difference schemes take this into account by upwind 
differencing at supersonic points and central differencing at subsonic 
points. Details differ for various schemes on the means of switching 
from one type of differencing to the other when the flow changes from 
supersonic to subsonic. 

A novel scheme was recently introduced by Moretti (Ref. 32) who was 
inspired by the physical implications of the method of characteristics. 
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The scheme is a second order accurate method capable of tracking weak 
shocks and is currently being improved, with initial testing being 
conducted mainly on unsteady problems in two space dimensions. It is 
appropriately called the X-scheme and employs shock fitting to specify 
the entropy jump and to correct the shock speed for shocks that are 
initially predicted as being isentropic. Other investigators have made 
contributions to this technique, in particular Zannetti (Ref. 33) who 
introduced the use of generalized Riemann invariants. Still, while the 
method offers hope of accurately tracking and computing shocks, its 
forthwith use for three-dimensional flows has yet to be demonstrated. 

The difficulties in using shock fitting procedures for complicated flows 
encouraged efforts for alternative methods. 

An alternative method which models the dissipative nature of shocks 
was introduced by von Neumann and Richtmyer (Ref. 34). Their approach 
was to add terms, mocking the physical effects of viscosity, to the one- 
dimensional unsteady inviscid flow equations. This was done in a 
discriminatory fashion by prescribing the additional terms as nonlinear 
functions of the dependent variables which have the following three 
features: to act in a strong manner at a shock location; to have a weak 

effect away from the shock; and to remain compatible with the Rankine- 
Hugoniot relations for shocks sufficiently thin when compared to the 
characteristic flow dimensions. The overall effect of this approach was 
to smear any shock smoothly over an interval of distance rather than 
produce a sharp discontinuity. The thickness of a typical smeared shock 
is on the order of the grid size interval which theoretically can be as 
fine as desired, though practical considerations limit the spacing 
refinement. The thermodynamic variables therefore acquire a continuous 
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but rapid variation in place of a sharp jump across a shock. Since 
these terms are arrived at based only on considerations of physical 
reasonableness rather than exact laws, a certain amount of arbitrariness 
is introduced; this would likely be compounded in generalization to 
higher dimensions. It should also be emphasized that the modifications 
introduced by including these dissipative terms is purely a mathematical 
artifice and the success of the technique relies on properly specifying 
the amount and distribution of artificial viscosity. The inclusion of 
physically reasonable, yet artificial, terms to provide treatment of a 
shock discontinuity illustrates a prevalent circumstance appearing 
variously in computational applications; namely, the acceptance of some 
loss in one property for a gain in another. Here this occurs with a 
highly simplified treatment of shocks being obtained at the expense of 
introducing some numerical error into the solution by using somewhat 
arbitrarily defined dissipative terms. 

The introduction of artificial viscosity preceded the advent of the 
high speed computer but was recognized, along with other developments in 
numerical methods, as being most useful with automatic computing 
capability. The rapid growth in computer capabilities signaled a trend 
toward application of highly simplified procedures suitable for 
repetitive machine calculation. The stage was set for numerical 
calculation of nonlinear flows with shocks. The amount of numerical 
work dealing with transonic flows is overwhelming and precludes any 
complete review here. A number of excellent surveys of developments in 
transonic flows have been made (Refs. 35 to 39. Therefore, only a brief 
review will be given to highlight certain key aspects. As in the case 
of analytical treatment of propeller theory, it is useful and convenient 
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to sketch the numerical approach by focusing on tne numerical 
developments associated with wing applications. Most of the fundamental 
advances in the numerical treatment of external transonic flows have 
come about largely due to procedures devised to predict steady and 
unsteady flows over airfoils. These flows possess the characteristics 
listed above and yet offer examples simple enough to be treated directly 
with various techniques. This quality affords the opportunity to 
compare the numerical results and lend credibility to them; also, 
experimental data or exact analytical solutions for restricted cases can 
often be obtained to qualify or corroborate the numerical results. 

An appropriate beginning for the discussion of transonic flows 
over airfoils is the work of Magnus and Yoshihara (Ref. 40). This work 
represents the first direct finite-difference procedure for transonic 
flow with embedded shock waves. The unsteady Euler equations, expressed 
in conservation form, were used to solve two-dimensional, steady, 
supercritical flows about lifting airfoils. These unsteady equations 
were numerically approximated using a modified Lax-Wendroff finite- 
difference technique. Steady solutions were obtained by allowing the 
unsteady solutions to converge to a time-dependent state. The advantage 
in treating a steady flow as a limit of an unsteady one is that the 
steady equations are a mixed elliptic-hyperbolic system whereas the 
unsteady equations are hyperbolic. The gain realized from this 
uniformity of equation type more than offsets the increase in complexity 
caused by adding time as another independent variable. This approach 
especially simplified the problem since the locations where the steady 
flow equations switch type are unknown and must be found as part of the 
solution process. Use of the unsteady formulation eliminates the need 
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to devise distinct numerical schemes particular to the subsonic and 
supersonic regions. Of considerable importance was their method of 
introducing numerical viscosity. They introduced artificial viscosity 
in two separate ways. One of the ways diffusion was introduced arises 
from the formulation of the second order differencing scheme. Lower 
order terms in the truncation error are retained which behave similarly 
to real viscosity by reducing all gradients in the solution. The higher 
gradient locations naturally were affected more than the lower ones. 

This serves the same purpose as the artificial viscosity introduced by 
von Neumann and Richtmyer in that it automatically captures a shock as a 
high gradient region joining regions of strongly dissimilar flows. A 
second way of introducing damping terms was brought about by using 
spatially weighted averages of the dependent variables in the leading 
terms of the Taylor series expansions expressing the updating of these 
variables. The numerical viscosity can be reduced in each case by 
refining the grid since the the damping terms are in fact truncation 
error terms. A trade-off is reached between stability and accuracy. 

The solutions obtained were quite time consuming, in good part due to 
the absence of a suitable grid, but served as benchmarks for more 
approximate calculations to follow. 

Murman and Cole approached the transonic flow problem along the 
lines of small perturbation analysis. They first solved a simplified 
form of the two-dimensional steady transonic equation 



for the perturbation potential <p with what was later realized as not 
being a fully conservative scheme (Ref. 41), and subsequently with a 
nearly fully conservative scheme (Ref. 42) in which the correct shock 
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jump condition was obeyed, at least in practice as evidenced by the 
results. The perturbation potential is the velocity potential minus the 
free stream potential. Here, and in what follows, the x coordinate is 
in the free-stream direction; the y direction is normal to the free- 
stream flow. The use of a perturbation technique limits these efforts 
to thin airfoils but greatly simplifies the boundary condition treatment 
at the airfoil surface. They used an iterative line relaxation 
procedure to sweep downstream along successive grid lines aligned 
transverse to the flow. Over-relaxation was used in the subsonic 
regions of the flow. The use of the successive line over-relaxation 
(SLOR) procedure improved the computational efficiency by about an order 
of magnitude over that of Magnus and Yoshihara. 

The use of such a straightforward procedure for the calculation of 
mixed flow was possible due to their type-dependent differencing scheme. 
This was their main contribution and is carried out as follows. During 
the course of the calculation, the velocity at each grid point along the 
relaxation line under consideration is compared with the speed of sound. 
If the velocity is supersonic, that point is implicitly backward 
(upwind) differenced in the streamwise direction. If the velocity is 
subsonic, central differencing is used. This type-dependent 
differencing respects the domain of dependence and furthermore, proved 
simple to implement. An important development was included in the fully 
conservative scheme which significantly improved the shock calculations. 
A special difference operator was introduced which was used at grid 
points where the flow decelerates through the speed of sound (ruling out 
any supersoni c-to-supersoni c shocks). The difference operator thus 
operates at shock points and is called a shock-point operator. It 
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should also be mentioned that a switching operator (Ref. 43) is used at 
sonic points, where the flow accelerates through the speed of sound. 

The sonic operator poses no special problem since the flow is continuous 
at sonic points. However, it should be emphasized that four operators 
exist as fellows: a central operator at subsonic points; an upwind 

operator at supersonic points; a sonic operator; and a shock-point 
operator. The shock-point operator effectively switches the difference 
scheme from the implicit backward hyperbolic operator upstream of the 
shock to the centrally differenced elliptic operator past the shock in 
such a manner that the combination of all three operators is fully 
conservative . 

This switching procedure can be looked at from a flux viewpoint by 
imagining the grid points to be enclosed by corresponding cells packed 
tightly. Use of the switching operator balances the flux internally so 
that for any contiguous group of grid cells, the net flux of a quantity 
at their external border is zero. Since the initial type-dependent 
scheme lacked the shock-point operator, flux was not necessarily 
conserved at the internal shock boundaries and consequently the scheme 
was not fully conservative even when the derivatives were otherwise 
placed in conservation form. 

Several codes have originated based on the Murman and Cole method. 
The first one developed was by Krupp and Murman (Ref. 44) which used the 
nonconservative form of the difference equations. Murman, Bailey, and 
Johnson (Ref. 45) later developed a code (TSFOIL) using the fully 
conservative form of the difference equation given by 

(' - H -K - *4 * fyy - 0 


(3.25) 
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or, in a general form, with the physical variables scaled, written as 


( K *« - H - 1 ? x) x * ■ °- 


(3.26) 


The relations between the physical and scaled variables are given as 
fol lows : 


where 


--2/3 M n 
9 = 6 M cp 


y = s 1/3 m^ 


K = 



M 2m S 2/3 


(3.27) 


The variable S is the airfoil thickness ratio, is the free- 
stream Mach number, and H is the ratio of specific heats. The 
constants d, m, and n are the similarity parameters chosen to give 
the desired form of the small disturbance equation. When d = m = n = 0, 
the Cole form (Ref. 46) is obtained. The Sprieter form (Ref. 47) is 
given by defining d = 2 and m = n = 2/3. Finally, the Krupp form 
(Ref. 48) is obtained by selecting d = 7/4, m = 1/2, and n = 3/4. The 
Cole choices for the transonic similarity parameters are based on 
mathematical simplicity. The Spreiter and Krupp choices are made to 
give accurate approximations over a wide range of free-stream velocity 
and thickness ratio. The conservative form of the equations produce 
shocks with nearly the correct strength and speed. Nonconservative 
formulations exhibit the fortuitous feature of having the shock position 
displaced upstream accompanied by a smaller pressure drop across the 
shock, more in agreement with experimental viscous measurements. Such 
features make the nonconservative results more appealing to some airfoil 
designers . 
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A few loopholes present in the Murman-Cole shock operator were 
closed with the simple modifications made by Engquist and Osher 
(Ref. 49). They not only analyzed the two-dimensional small disturbance 
equation, Eq. (3.26), but solved the time-dependent form: 

- ^-f- 1 ^ - 2$ xt « 0. (3.28) 


Convergence to the steady state was reached after a sufficient number of 
iterations . 

Artificial viscosity was not added by design in the numerical 
formulations of the type dependent schemes mentioned above. However, 
enough is generated due to the upwinding operators that shocks are 
captured over 3 or 4 mesh widths using the Murman-Cole scheme and over 
1 or 2 mesh widths using the Engqui st-Osher scheme. The use of type 
dependent differencing destroys the symmetry that would otherwise occur 
in the equations and which would admit expansion shocks as well as 
compression shocks. 

The SLOR procedure used by Murman and Cole was replaced by 
approximate factorization (AF) procedures by Ballhaus and Steger 
(Ref. 50) in their investigation of the low frequency form for the 
transonic small disturbance equation which is given in the form 


2 


03 0 . 

U 


00 


M 2 <p 

no' 


xt 


V <p 
c T xx 


+ <P' 


yy 


(3.29) 


for an airfoil of chord ft. undergoing periodic motion of frequency co. 

2 m 

The free-stream velocity is and V c = 1 - - (H + 1 )M 00 < P X - The 

flow is locally subsonic or supersonic when V c is positive or negative, 
respectively. The ratio coft/U^ determines the time scale and is called 
the reduced frequency. It relates the time a mean fluid particle takes 
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to traverse the length of a chord to the time required for one blade 
oscillation. The appearance of the reduced frequency is attributed to 
the fact that in Eq. (3.29) the variables 9, x, y, and t are 
additionally scaled by Q-U^, a, a, and 1/to, respectively. This low 
frequency form of the equation is valid for reduced frequencies much 
less than one. For most transonic flows the low frequency equation is 
adequate since the low frequencies associated with the slowly moving 
upstream waves dominate the solution near the airfoil. The higher 
frequency waves travel quickly downstream and do not have a chance to 
build into strong waves. Fortunately, the low frequency oscillations 
are adequately analyzed with coarser time steps than those for high 
frequencies. For this reason Ballhaus and Steger examined implicit and 
semi-implicit numerical schemes that would allow them to remove the 
time-step limitations of explicit schemes so that larger time steps 
could then be taken which were limited only by accuracy considerations. 
The AF procedures were introduced to reduce the computational work 
required by SLOR techniques. The AF procedures split the difference 
equations into simple factors which are easily invertible. The factors 
are usually solvable by tridiagonal or quadridiagonal algorithms such 
as the Thomas algorithm. Among the schemes studied was the ADI method 
developed by Peaceman and Rachford for the two-dimensional Laplace 
equation and extended to multi-dimensions by Douglas and Gunn (Refs. 51 
to 52). The ADI scheme was shown to be an AF type scheme. The 
stability of the schemes was compared and their shock capturing features 
were examined. 

The ADI technique developed by Douglas and Gunn was used in the 
numerical solution of the unsteady, three-dimensional flow about 
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helicopter rctors (Refs. 1 to 5) using the small disturbance approach. 

A numerical code was written and used to analyze the flow in the 
transonic tip region of the rotor blades. The strong similarity in the 
form of the helicopter equation to that of the propeller problem led to 
the adoption of the code, with some modifications, for the present 
effort. 



IV. GOVERNING EQUATION 

This chapter presents the potential equation and the development 
of an approximation to it which is valid for small disturbance flow 
about thin, lightly loaded propeller blades. First, the general 
potential equation in an inertial coordinate system will be developed 
for isentropic flow of a perfect gas with no other approximations made. 
Then this equation will be transformed to streamwise coordinates. The 
potential equation in rotating coordinates will then be given as 
developed in (Ref. 5) for a coordinate system attached to a helicopter 
blade. Following this, the equivalent tensor equation valid for any 
coordinate system will be presented for an accelerated system. Finally, 
the small perturbation equation for generalized helical coordinates will 
be derived. 

4.1 Potential Equation in an Inertial Coordinate System 

Consider a flow with uniform upstream velocity directed 
along the X-axis where X, Y, and Z are Cartesian coordinates with 
corresponding velocity components u, v, and w. The continuity 
equation can be written in the expanded form as 

Pj. + UP^ + Vpy + Wp^ + p(Uy + Vy + w^) = 0. (4.1) 

Bernoulli's equation for unsteady flow can be expressed in a form 
giving the speed of sound c as a variation from its upstream speed 
c as 

00 

c 2 = c 2 - (24> t + 4> x + 4>y + 4> z ) (4.2) 
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where the scalar <1> is the velocity potential defined through the 
following components: 


u 


3$ . 

' 9X ’ 




(4.3) 


Thus, the fluid velocity is given by 


q = V4>. 


(4.4) 


For isentropic flow, the density and speed of sound are connected 
by the relation 


= 


p 


CD 


feP” 


(4.5) 


This relation can be used to eliminate the density in Eq. (4.1). Using 

Eq. (4.5), performing the differentiation in the first four terms in 

2 

Eq. (4.1), and multiplying the result by c /p yields 


j_ r (c 2) 

+ 1 [ t 


+ u(c ) x + v(c ) Y + w(c ) z + c (u x + v y + w z ) = 0. (4.6) 


Bernoulli's equation, written for an unsteady flow in Eq. (4.2), can be 

2 

used to introduce the velocity potential into Eq. (4.6); replace c 
in the equation using Eq. (4.2), and carrying out the indicated 
differentiation yields 


*tt + “"Vt * v *Yt * “*Zt * U<<t xt * * V ^XY + "V 


+ V(4>y t + U4> X y + v4>yy + W<I>y 7 ) + W^j. + + V*^ + W^) 

- c 2 (4> xx + 4>yy + 4> zz ) = 0. (4.7) 


Rearranging terms yields a more symmetrical equation, viz., 

2 2 2 2 

<t>^ + 2u 4> x ^. + 2v4>y t + 2w 4> z ^ + (u - c )4> xx + (v - c')4>yy 

9 2 

+ (w - c )4> zz + 2uv4> xy + 2uw$ xz + 2vw4>y Z = 0. 


(4.8) 
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Various approximations can be made to Eq. (4.8) which itself is an 

2 

exact equation. For example, dividing it by c and allowing c to 

2 

approach infinity reduces the equation to V 4> = 0 for an 
incompressible flow. Similarly, other approximations can be made by 
removing cubic and even quadratic terms. As another example, in the 
case of acoustic theory, a linear equation results for the subsonic and 
supersonic regimes by neglecting certain terms and letting the 
quantities c and u assume their free-stream values c^ and U^; 
whereas a separate equation is obtained for the transonic regime since 
local effects are more dominant, requiring retention of the local values 
of c and u. Most approximate equations derived from Eq. (4.8) are 
produced for the purposes of linearization. However, the goal here is 
to arrive at an approximation that is consistent with this equation in 
the limit of small disturbances, yet valid for nonlinear transonic 
flows . 

4.2 Potential Equation in Streamwise Coordinates 
Before beginning our discussion of the disturbance potential for 
the propeller problem, it is worthwhile to present Eq. (4.8) using 
streamwise coordinates. First, the magnitude of the total velocity 
vector q is given as 

q - (u 2 ♦ v 2 * w 2 ) 1 ' 2 - <4 - 9> 


For any directional coordinate s, the directional derivative along s 
is simply 


(4.10) 
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By referring to Fig. 4.1, It can be seen that dX/ds, dY/ds, and dZ/ds 
are the direction cosines which equal u/q, v/q, and w/q, respectively, 
If s Isa streamline. The second directional derivative Is defined as 

9 2 $ 1 3 2 4> 2 9 2 4> 2 A ? 

as 2 ss q 2 V ax 2 ax 2 az 2 3X3V 


+ 2uw 


8 2 4 > 

9X9Z 


2vw 


9 2 4> 

9Y3Z 


(4.11) 


Furthermore , 


8 2 4, 1 ( u . v + w ak. 

9 S 3 1 q l U 3X3t + V 3Y3t + W 3Z9t 


(4.12) 


follows directly from Eq. (4.3) when the velocity ratios are Inserted as 
the direction cosines. Substitution of the streamwlse derivatives given 
by Eqs . (4.11) and (4.12) Into the potential Eq. (4.8) Is easily 
accomplished by grouping the appropriate terms In the latter. This 
leads to a streamwlse formulation of the potential equation given by 


4> tt + 2q<t> st = (c 2 - q 2 )4> ss + c 2 (V 2 4> - $„) 


(4.13) 


st ■ "* ^ '’ss ss 

where V 2 = 9 2 /3X 2 + 3 2 /3Y 2 + 9 2 /3Z 2 Is the Cartesian Laplaclan 
operator. This takes on a particularly simple form for two-dimensional 
steady flow, 1 .e . , 


(c 2 - q 2 )* ss ♦ c 2 * nn - 0. 


(4.14) 


where 4> ss reduces to 


. 1 / 2 9 2 4> , 11W 9 2 4> 2 

“ q 2 V ax 2 8X3Y ax 2 


(4.15) 


and 4> represents the second directional derivative normal to the 
nn 

streamwise direction such that the positive direction of n is 
cons 1 stent with Fig. 4.1: 




/ 


X 

FIGURE 4J. - DIRECTION COSINES EQUAL u/q, v/q, AND w/q 
WHEN s IS DIRECTED ALONG q\ 
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- 2uv 


3 2 <l> 

3X3Y 


+ u 


3 2 <E 

SY" 


(4.16) 


The use of a coordinate system aligned with the streamwise and 
normal directions has obvious advantages. The potential equation is 
considerably simplified when expressed in these coordinates. The 
troublesome cross-derivative terms which generally complicate the 
numerical solution and usually adversely effect the convergence rate are 
removed. Even when streamwise and normal coordinates cannot be used 
explicitly, they can be expressed locally in terms of the actual 
coordinates as in Eqs. (4.15) and (4.16) for two-dimensional flow. One 
marked advantage of doing this is that derivatives can be treated 
differently in the streamwise and normal directions. In fact, it was 
for this very reason that Jameson (Ref. 53) introduced his rotated- 
differences procedure. The use of a streamwise coordinate system allows 
type dependent differencing schemes, as discussed in the last chapter, 
to be employed in a straightforward manner. The streamwise direction 
can be type-dependent differenced and the normal direction can be 
centrally differenced. The use of rotated differencing allows 
derivatives to be constructed in the streamwise and normal directions 
from the Cartesian coordinates. Failure to use a rotated differencing 
scheme without one coordinate being nearly aligned with the flow 
direction can lead to instability when supersonic regions are present. 
The likelihood of instability results from either central differences 
having a component in the streamwise direction or from upwind 
differences which may not contain the correct domain of dependence, 
leading to negative artificial viscosities. This points out a major 
advantage in using helical coordinates for propeller flow. The 
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undisturbed flow direction can be aligned with a helical coordinate. 

This helical direction will be the streamwise direction for which type- 
dependent differencing can be used. 

4.3 Potential Equation in Noninertial System 

The potential equation, Eq. (4.8) can be expressed in vector form 
(Ref. 54) as 

$ tt + ft (v4>)2 + * v(v<1))2 ■<£-<*- 1J [*t + 1 <V<J>)2 

The potential equation can also be expressed in an xyz coordinate 
system which is both translating with an arbitrary constant velocity V 
away from the inertial frame, and rotating with an arbitrary constant 
rotational velocity fi x *. The angular velocity vector Q (with 
magnitude Q) has the direction of the axis through the origin about 
which xyz are rotating, and % is the position vector in this 
system. The relationship between the two coordinate systems is shown 
in Fig. 4.2. Although both V and Q are general vectors, they will 
soon be restricted to a common axial direction in the following 
development. The rotating system is, of course, not an inertial system. 
Although such a system introduces extra terms, it offers advantages over 
fixed coordinates. In particular, for a propelle*" rotating about an 
axis aligned with its flight direction, the flow appears steady 
(ignoring unsteady effects such as flutter) to an observer in a 
coordinate system rotating with the propeller. 

A potential equation can also be established in the noninertial 
frame for a flow which is irrotational in the inertial frame. This is 
initiated by defining a perturbation potential cp which separates the 
contribution of the free-stream from the total potential. The 




(4.17) 
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FIGURE A. 2. - RELATIONSHIP BETWEEN INERTIAL COORDI- 
NATE SYSTEM XYZ AND NON INERTIAL COORDINATE SYSTEM 
xyz WHICH IS ROTATING ABOUT AXIS PASSING THROUGH 
ITS ORIGIN. 
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freestream contains both uniform flow and bulk rotational flow 
contributions, consequently 


q' = Vcp = q - (V + Q x J) (4.18) 

represents the irrotational perturbation velocity. The gradient 
operator is invariant with respect to choice of Cartesian coordinate 
systems. The components of the perturbation velocity q' are defined 
as 


u « = & 
u 8x 


v 1 = ^ 

3y 


w 1 - 9* 

W ' 8z 


(4.19) 


The transformation of the potential equation has been carried out by 
Isom (Ref. 5) for a coordinate system translating and rotating with 
constant, but otherwise arbitrary, V and 0. The potential equation 
in the noninertial system is given as 


<p.|.£ + a • V(V«p • a) + 2a • V<p^. - Qxa • Vcp - • 7cp + 27<p • 7<p^ 

+ a • 7(7<p) 2 = jc 2 - a - 1 )^9 t + 

where a is the negative of the vector sum of the translational 
velocity and the rotational velocity as seen by ar. observer in the 
rotating frame. Thus, for such an observer, a can be simply written as 


a • 7<p + l(7<p) 


% 


7 2 <p 


(4.20) 


a = -(V + Q x 4). (4.21) 

It should be pointed out that in the transformation from the 
potential equation for an inertial reference frame, Eq. M.'/y. */, t/c. 
potential equation for a noninertial frame, Eq. (4.20), the cubic term, 
on the left-hand side of Eq. (4.17), was dropped. Otherwise, Eq. (4.20) 
is an exact equation. Introduction of the perturbation potential does 
not, in itself, imply any approximations. 
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4.4 Potential Equation in Noninertial Curvilinear Coordinate System 
Application of Eq. (4.20) to the propeller problem can be carried 
out by further simplification in much the same manner as was done in 
deriving the small disturbance equations for the helicopter problem. 
However, since the helical coordinate systems which will be used for the 
propeller problem are nonorthogonal , the potential equation, Eq. (4.20) 
is more tractable if expressed in invariant tensor form. For a 
translational velocity V which is independent of time, it has the form 


' - L *-- J * , 3 k <v k - v,> % 

at L J J ay 1 


&♦»’ ^V 2 *’ ^ 

at 2 ay \ ay 3 I ^ 

, i. (ji 3sl a* X K . | <£ - (t - 1) 
3t l 9 ay' ayjj ay’ V ^ */) 


a^ i 3 c£_ 
9t ay 1 


. I JJ 3je_ SjL. 

2 9 3y' ay 3 


Ij a 2 <p j_ x (i/gg 53 ) ^ 
9 ayV V’ay’ V «»' 


(4.22) 


In this equation a 1 and a. represent the contravariant and covariant 

components, respectively, of the general velocity vector a defined by 

Eq. (4.21). The quantities g 9 are the contravariant components of 

the metric tensor for the transformation between the curvilinear 

coordinates y^ and our orthogonal Cartesian coordinates. The quantity 

g is the determinant of the corresponding covariant components g^ 

i i k 

which do not appear explicitly in the equation. The symbol e is 
the permutation symbol which equals zero for repeating values of i.j.k, 

I 

i unity for cyclic (even) permutations of 1,2,3, and negative unity 

otherwise. 

4.5 Potential Equation in Helical Coordinates 
Equation (4.22), expressed in terms of a rotating coordinate 

system, can be simplified considerably by choosing one of the 

i 


I 
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curvilinear coordinates to be in the direction of the vector a. From 
Eq. (4.21), this is seen to be the direction opposite the vector sum of 
the translational velocity and rotational velocity which we now take to 
be orthogonal such that V is aligned with the axis of rotation as 
shown in Fig. 4.3. With this arrangement the vector "a is in the 
direction of the free-stream velocity vector as it appears to an 
observer whose frame of reference rotates with the blade. In this 
instance the magnitude of aT will be designated by U such that it may 
be written in matrix form as 


Also, 



U 2 = (fir) 2 + V 2 , 


(4.23) 


(4.24) 


where r is the radial distance as measured from the axis of rotation. 
Equation (4.22) then simplifies to 


, u2 , 2U -A. 
8t Cay 1 ) ayat 


1 q 2 v 2 . 9_ / Ij 9f£L 

2 y 2 8t l y i i 

ay V ay ay J . 
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ik a<e_ 8 c£_ 
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n x || 9$_ . 1 Jj 3 sl I 
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V 5 3y 

4.6 Approximate Potential Equation in Scaled Helical Coordinates 
Equation (4.25) is exact in that it is equivalent to Eq. (4.20), 


(4.25) 


but expressed in a blade-fixed reference system. However, this equation 
is much too complicated to be solved efficiently. Therefore, a 
systematic simplification is necessary to arrive at an approximation 
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which retains the nonlinear features of the problem and which is valid 
for small disturbances about the mean flow. Starting from Eq. (4.20), 
such a procedure was carried out by Isom (Ref. 5) using Cartesian 
coordinates for flow about the tip region of a helicopter blade. 
Following the development of the approximate equation for the helicopter 
problem, a similar process is applied to Eq. (4.22) for the case of a 
propeller blade using helical coordinates. The use of helical 
coordinates introduces the metric tensor into the approximate equation 
and allows representation of flow curvature. 

Derivation of the approximate equations for flow about a propeller 
is based on consideration of the following parameters: the axial 

free-stream Mach number M ; the thickness ratio 5; the ratio of 
the chord 2. to the blade-tip radius R, which is the inverse of the 
aspect ratio, c; and the advance ratio X. For true unsteady problems 
the reduced frequency would also enter. However, since we are examining 
only steady solutions, the reduced frequency need not be given further 
consideration. Also, for cascade solutions, a parameter representing 
the blockage of the flow, say solidity, would normally enter the 
problem; here the blades will be considered far enough apart so that 
variations in solidity will not be important. For an advanced 
turboprop, typical values of the relevant parameters are = 0.8, 

8 = .02, e = 0.3, and X = 1 . 

Following Cole's choices of the similarity parameters, the lateral 

direction y , which lies essentially normal to the blade, is scaled 

1 / 3 

by dividing it by the value 8 . This could be done also for the 

2 

radial-like or spanwise direction y , since the physical justification 
for scaling is to transform the lateral dimensions to account for the 
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weak diminution of disturbances in these directions in comparison to the 
streamwise direction. However, Isom chose to scale the spanwise 
direction differently by introducing e as defined previously. This is 

accomplished by normalizing the streamwise direction y 1 along with the 

^ 2 
already scaled y direction by d, while normalizing y by R. 

Since by these two methods of scaling, the y 1 result in the same final 

approximation, either can be used. The latter choice is made here and 

the directions are scaled such that the original coordinates transform 

to the dimensionless ones, as given by 

y'.ty', y 2 = Ry 2 . y 3 .-^? 3 (4.26) 

6 

where the tilde denotes the dimensionless coordinates. In addition, 
time and the disturbance potential are nondimensional ized as follows: 

t = ^ , <p = S 2/ W (4.27) 

The use of the tilde is for clarification only. When the switch to the 
scaled, nondimensional variables is actually performed, the tilde will 
be removed. 

It is now assumed that in terms of the scaled variables, the 
following are two reasonable approximations: (1) all second derivatives 

of the potential are of first order in magnitude; and (2) all first 
derivatives are second order in magnitude. The first assumption is 
based on accumulated experience and the second is based on the 
fundamental limitation that the deviation in velocity from the free- 
stream velocity be small. For the limiting case of incompressible flow, 
the maximum variation from free-stream conditions occurs at the boundary 
of the solid surfaces; this can be expected to extend reasonably well to 
compressible flow. The perturbation velocities along a blade surface 
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can be considered small for sufficiently thin blades with small slopes 
and thus all perturbation velocities can be considered small. Then, In 
terms of the original unsealed variables, using Eq. (4.26) and the 
assumption of equal magnitude of scaled second derivatives, the 
following ratios among the second derivatives, 

o ? ? 2 2 2 

3 <p . 3 y . 8 <p . 8 <p . 8 <p . 8 <p 

(ay 1 ) 2 ' (ay 2 ) 2 ' (ay 3 ) 2 ’ 3 »V ' a »' ; ^ ’ 8 *V 
and, similarly, the following ratios among the products of the first 
derivatives 


8<P \ 2 . ( 8<p \ 2 . / 8<p . 8<p 89 . 8<p 8<p . 8<p 8<p 


1 2 
8y 8y 


1 3 * 2 3 

8y 8y 8y 8y 


V/ \8y 2 / \8y 3 

are each, in the given sequence, estimated to compare as 
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1/3 


1/3 
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The ratio relations above can be multiplied by 2. and written 

. 2 . .2/3 . . ,1/3 . .1 /3 

1 : c : 8 : e : 8 : eS 

From the values given earlier for a typical advanced turboprop blade, 
e 2, 6^3, and are each small compared to unity. Furthermore, 

the values of the off-diagonal elements of the contravariant metric 
tensor components g 1 '-' are small compared to the diagonal elements 
for most locations, especially near the blade surfaces where 
gradients in the solution are largest. These are the main criteria used 
to reduce Eq. (4.25) to a more amenable form. 

The expression g 1 ^ 3<p/3x' 8cp/8x^, representing (V<p) , occurs in 
three terms of Eq. (4.25) and is approximated as 

gi j 3^4 3^ * g 1 1 / 

3x^ 3x^ \ 3x 1 / 


(4.28) 
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This is justified since each of the discarded terms is judged to be 
small by virtue of the corresponding ratio relations and by the fact 
that the off-diagonal g 1 ^ are small. Of these three occurrences, 
only its occurrence in the last term on the left-hand side of Eq. (4.25) 
will be retained. This term may be written as 


g n/^\ 2 l,u^/ 2g n 

L l »»V J ay' \ ( 3y ') ay’ ) 


(4.29) 


The variation of g 11 with respect to y 1 is small enough that the 
second term on the right-hand side of Eq. (4.29) can be neglected. This 
results in the following approximation 


u 3 g n/vy sUg ii ji/asV 

3y V 3y / 3y \3y / 


(4.30) 


The Laplacian operator occurring in the final brackets on the 
right-hand side of Eq. (4.25) remains to be simplified. This is done by 
neglecting the last set of terms which contain the first derivatives of 
the potential. The Laplacian factor is thus given the approximate form 


s’ 3 9 ij ) v* 

3y' 3y^ T y 3y 3Y J 


ij _jr 


3y ] 3y J 


(4.31) 


Equations (4.30) and (4.31) are inserted into Eq. (4.25) and the fourth 
term containing a first derivative of the potential on the left-hand 
side of the equation is omitted. This reduces Eq. (4.25) to 
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(4.32) 
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Finally, on the right-hand side, the unsteady term is dropped and only 
the largest of the nonlinear terms is retained which contains the factor 
g^ 3 2 cp/3y*3y^ from the Laplacian. All other nonlinear terms are small 
based on the equivalence of the ratio relations. The equation for the 
perturbation potential is then 


^ ♦ u 2 . 2U ♦ 2Ug' ' - 1 

3t (ay 1 ) 3y 3t 3y Cay’) 


= -a - nug 11 a<p ? + c 2 

3y (ay 1 ) 2 


11 


+ 9 


22 3 cp 


(ay 1 ) (ay 2 )' 


+ g 


+ 2^g 

By defining 
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(4.33) 



(4.34) 


and combining similar terms yields, after dividing by 
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(4.35) 


3y'3y" 3y ay* 

This is the final form of the perturbation potential equation in 


unsealed variables. This equation will now be transformed by the 
scaling relations given in Eqs. (4.26) and (4.27). Substitution of 


these scaling relations into Eq. (4.35) gives 
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( 4 . 36 ) 


2 2/3 

where the tilde is not shown. Multiplying this equation by i /$ 
and introducing a constant Mach number, characteristic of the rotational 
speed of the propeller tip as defined by 


M - 
T c^ * 


( 4 . 37 ) 


results in 
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( 4 . 38 ) 


For low frequency or steady problems, it is admissible to delete 
the second derivative with respect to time. This leaves the final form 


of the governing equation in terms of the scaled variables as 
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S 8y 3y* 

Equation (4.39) is the small disturbance equation which is solved 
numerically by the ADI Douglas-Gunn algorithm presented in Chapter VII, 



V. BOUNDARY CONDITIONS 


5.1 Wall Tangency Condition for Blade Surfaces 
The small disturbance boundary condition to be applied at the blade 
surfaces will now be established. This necessitates returning briefly 
to the unsealed variables. 

In a blade-fixed coordinate system, if the equation for the surface 
of a blade moving in a time dependent manner is given by 

F(y 1 , y 2 , y 3 ,t) = 0, (5.1) 

then the vanishing of the fluid velocity component normal to the surface 
brings 

FCy 1 , y 2 , y 3 , t) = || + (a + V<p) • VF « 0 (5.2) 

where 

F = F u (y\ y 2 , y 3 , t) = y 3 - h^y 1 , y 2 , t) « 0 (5.3a) 

on the upper (suction) blade surface from the leading edge (L.E. ) to the 
trai 1 ing edge (T. E . ) , and 

F = F o (y\ y 2 , y 3 , t) = y 3 - h^y 1 , y 2 , t) = 0 (5.3b) 

on the lower (pressure) blade surface. Here, h is the profile 
parameter. This follows the convention used for wing surfaces. The 

coordinates are aligned as shown in Fig. 5.1, where y^ is nearly 

3 2 
along the mean chord, y is nearly normal to the chord, and y 

lies in the spanwise direction. 

The general time-dependent boundary condition which allows the 
profile parameter to include pitching, bending, and twisting is 
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FIGURE 5.1. - ARRANGEMENT OF HELICAL 
COORDINATES WITH RESPECT TO PRO- 
PELLER CASCADE. 
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developed next. This time-dependence is included mainly for any future 
efforts which consider flutter. Proceeding with this general 
development, the two spatial terms in Eq. (5.2) are 

9F 


a • VF = U 


3y 


1 


(5.4) 


and 


V<p • VF = g 


ij 9£_ 9F_ 


9y ] 3y J 


(5.5) 


The latter expression can be written out using Eq. (5.3) as 

n OF ii a*_ ah 22 a^_ ah_ 33 a$_ 
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(5.6) 


The usual application of small disturbance boundary conditions requires 

that an approximation be made to the derivative of the potential in the 

direction normal to the mean chord. The use of Eq. (5.6) along with 

Eq. (5.4), when substituted into Eq. (5.2), with rearrangement to 
3 

separate a<p/3y , gives 

2*, /g33 13 ah 23 ah \ |h „ y 3h + g ll 3^ 3JL. , g 22 5!L 

9 y 3 \ 9 ay 3y 2 / 3t 3y ay' 3y ay 2 3 T 

12 / a<£_ 3h _ 3cp 3h \ 13 3ce_ _ 23 &£_ 

+ g I 1 n + o’ T I - 9 1 9 0 • 




ay' ay" ay 2 sy 1 / 9 y' 

Introduction of the scaling laws into the above equation provides 
a basis for simplification. The scaling of h that is consistent with 
the scaling relations introduced previously is 

h = Q-StKy 1 , y 2 , t) . 


(5.7) 


(5.8) 


91 


Substituting Eq. (5.8) along with the previous scaling relations into 
Eq. (5.7) and again discarding the tilde, yields 
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(5.9) 


Based on the typical values of e and 6, this is closely approximated 
by 
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(5.10) 


For blades which taper gradually toward the tip, the last term in 

Eq. (5.10) can be neglected with the result 

9<2 ]_ ( 8h M 9h 13 1 9<p \ 

' ' _33 e 3t + M t - 1 9 .1/3 «J 1 ' 


(5.11) 
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One of the key developments of this work is the generation of 

1 3 

coordinate systems possessing the characteristic that g *0 at 
the blade surfaces. For these coordinate systems, the unsteady boundary 
condition at the blade surface, in scaled variables, is 


9<p 

3 

9y 
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9h M 3h 


,35 \ c at * « T , y i 


(5.12) 


It is interesting to note that in dimensional Cartestian coordinates 
1 2 3 

(y = x, y = y, y = z), this expression is the familiar 


89 3h 11 9h 
9z ~ at 9x • 

Finally, the steady scaled boundary condition is simply 

89 1 M 8h 

8y 3 * g 33 M T 8y ] ' 


(5.13) 


(5.14) 
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This boundary condition is not actually applied at the airfoil surfaces 

but rather at two surfaces separated slightly such that they straddle 

3 

the mean chord surface which itself lies on, or near, the y = 0 
surface. Thus, at the upper surface 


8c£_ 

ay 3 
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g 33 H 


3h 
T 3y 
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(5.15a) 


and, at the lower surface 
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3y 
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3h 
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(5.15b) 


5.2 Flow Field Boundary Conditions 


Additional boundary conditions must be specified to fully define 
the problem. These include the boundary conditions upstream and 
downstream of the propeller, along with those for the inner and outer 
radial surfaces, and those far above and below the blade. Additionally, 
the Kutta condition at the trailing edge of the blade is required along 
with a treatment for the trailing vortex sheet. The handling of a 
cascade of blades is identical to that of an isolated blade aside from 
the specification of the conditions above and below the blade. Even in 
the case of cascade flow, the flow region about a single blade is 
sufficient provided the application of periodic conditions reduces the 
extent of the problem. As mentioned previously, shocks require no 
special boundary treatment as internal discontinuities because they are 
captured numerically as part of the solution process. 


The upstream boundary condtion on the perturbation potential is 

cp = 0. (5.16) 


93 


This merely states that the free-stream condition prevails sufficiently 
far upstream. 

In the downstream direction, the value of the perturbation 
potential is unknown, but at a large enough distance downstream 

= 0 ( 5 . 17 ) 

3y 

is used as the boundary condition. 

At the hub, or inner spanwise boundary, the condition 

^ = 0 ( 5 . 18 ) 

ay 

is used, and the outer spanwise boundary condition is simply 

<P = 0 . ( 5 . 19 ) 


o 

Although y is not strictly a radial direction, it closely 
approximates it for the coordinate systems employed here. 

The far field above and below the blade has two boundary condition 
options. For an isolated blade, the condition 


&= = 0 ( 5 . 20 ) 

ay 


is assumed to apply both far above and far below the blade. This 

3 

requires a coordinate system whereby the y direction is normal to 
the free-stream direction y^ at these boundaries. For multiple 
blades, cascade-type periodic boundary conditions are applied such that 


c P (y 1 , y 2 , t) 


= ?(y\ y 2 , t) 


y = 


upper 

surface 


lower 

surface 


( 5 . 21 ) 


is satisfied over the boundary surfaces which are chosen to lie at the 
midchannel. This, in return for straightforward application, requires 
a coordinate system which provides nodal coordinates having azimuthal 
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periodicity. Such periodic coordinate systems are discussed in a 
subsequent chapter. Only solutions using the cascade boundary 
conditions are presented in this work. 

It remains to discuss the boundary conditions in the vortex wake. 
Rollup of a wake is ignored and, the resulting vortex sheet is 
approximated as lying between the downstream extensions of the helical 
surfaces upon which the blade boundary conditions are imposed. The 
strength of the vortex sheet is assumed to be preserved as it convects 
downstream. In addition, the direction of the vortex vector is assumed 
to be "parallel" to the free-stream direction. Since the vortex sheet 
is a free surface, it cannot support a pressure difference. Along with 
these statements and the fact that the normal velocity across the sheet 
is zero, it follows that the tangential velocity jump occurring across 
the sheet is entirely in the spanwise direction. The boundary condition 
at the trailing edge can be written in terms of the scaled coordinates 
as 

[<p] T E = T(y 2 , t) (5.22) 

where T is the circulation around the propeller blade, and 

[cp] = ft) , - ft) 3 (5.23) 

w y =0+ y =0- 

is the jump in the perturbation potential across the wake. 


VI. HELICAL COORDINATES 


This chapter contains the development of special periodic helical 
coordinates suitable for propeller problems which include cascade 
effects. The coordinate transformation between this system and an 
orthogonal Cartesian system is specified so as to provide simple 
periodic boundaries, and to also provide orthogonal properties at the 
blades. However, before describing these new helical coordinates, a 
brief discussion will be given for a simpler helical coordinate system 
satisfactory for flow about an isolated blade, but not for cascade flow. 
This discussion points out similarities and distinctions between the two 
coordinate systems and will help to illustrate the features of the new 
coordinates. 

6.1 Helical Coordinates for an Isolated Blade 

Simple helical coordinates are useful in calculating the flow about 

isolated propeller blades. The transformation of these coordinates 

along with a listing of their metric tensors is given in Appendix B. 

These simple helical coordinates consist of two sets of helices, 

1 3 

providing two coordinate directions, y and y , with the remaining 

? 

coordinate direction y being essentially radial. A set of these 
coordinates is shown in Fig. 6.1. The helical curves are confined to 
spiral about a common axis on circular cylindrical surfaces; members of 
each set spiral clockwise about the axis when viewed along the axis in 
the direction corresponding to an increasing value of y\ The axial 
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i 3 

components of y and y are oppositely directed. The two sets 
of helices are mutually orthogonal. 

The main advantage of these coordinates is that the undisturbed 
flow direction can be aligned with one set of helices. These streamwise 
helices are constructed such that the helices of greater radii have 
smaller advance ratios; i.e., the smaller the radius the steeper the 
spiral, exactly characterizing the mean propeller flow. The set of 
helices normal to the streamwise helices behaves in the reverse way. 

Another distinct advantage is the orthogonality between the two 
sets of helices. This is an important feature because the small 
disturbance boundary condition applied at the blade surfaces is greatly 
simplified. Orthogonality is a desirable property because it generally 
results in more accurate numerical calculations. Thus, it is 
advantageous to have these coordinates orthogonal throughout the region 
of interest. 

2 

The coordinate direction y , which serves to measure the radial 
value is, in general, not a straight line. In fact, this coordinate 
direction has only one straight coordinate line: its coordinate axis, 

which was chosen to align with the leading edge of the blade. This is 
called the pitch change axis. With increasing distance from the pitch 
change axis, this coordinate direction departs to an increasing extent 
from a purely radial direction. Hence, the orthogonality of the 
coordinate system is reduced as the distance from the blade's leading 
edge is increased. Only at the leading edge is the coordinate system 
truly orthogonal. A notable charcter i sti c of the coordinate system is 
that, for high advance ratios, it tends toward an orthogonal Cartesian 
coordinate system. This can easily be seen by inspecting the metric 


t 
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tensor for this transformation which is provided in Appendix B. Since 

for all space g n - 1 and g 33 = 1, the values of y 1 and y 3 

for any given point correspond to the respective arc lengths along these 

2 

curves from the origin to that point. This is not the case for the y 
curves. As mentioned above, the y 2 coordinate is a measure of the 

2 

radial value from the axis to a given point; the arc length along y 
exceeds this value. However, as the advance ratio increases, g 22 
al so tends to unity. 

6.2 Periodic Helical Coordinate System for Cascades 

The periodic helical coordinates are closely related to the 

helical coordinates described above. The streamwise helices y 1 are 

unaltered so as to retain their alignment with the undisturbed flow. 

3 

However, to provide periodicity, the geometrical scaling of the y 
direction is modified so that this coordinate is no longer a direct 
measure of its arc length. Rather, the y direction is specified 
to scale similar to the angular coordinate of the familiar circular 

3 

cylindrical coordinate system; i.e., a change in the value of the y 

coordinate brings a change in arc length which is proportional to the 

radius. However, this change in scaling the y coordinate is not 

3 

sufficient to give circumferential periodicity. The y coordinate 
curves must be changed so that they are no longer helices. 

To provide circumferential periodicty, the y coordinate must 
be constructed such that, in a periodic fashion, it repeats its axial 
locations. The simplest choice is to make y the circumferential 
direction, forcing it to be independent of the axial value. The 
drawback in doing this is that then it would no longer be normal to the 
streamwise direction. Orthogonality in these two directions is an 
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important quality, especially near the surfaces of the blades. To 

3 

provide orthogonality at the mean blade locations, the y direction 
must assume the direction of the y coordinate of the helical system 
used for an Isolated blade, at least near the blades. Away from the 
blades, the y coordinate should reverse its axial direction so as 
to bend back and regain its original axial station. This must be 
repeated in a periodic fashion to conform to the locations of the 

blades. Thus, only coordinate curves in the streamwise direction y^ 

2 

are helices for this system. The y direction remains a coordinate 

which measures the radial value and since it is formed by the 

1 3 

intersection of surfaces of constant y and y , it will generally 
not be a straight line. 

A set of helical coordinates which has these properties is given by 
the following transformation to the Cartesian coordinates x 1 . 

= y 2 si n t? (6.1) 

x 2 = y 2 cos i? (6.2) 

x 3 = - q y 1 + A(y 2 )B(y 3 ) (6.3) 

where d is the circumferential angle measured as shown in Fig. 6.2, 
and which is related to the helical coordinates by 


d = U yl + R~ • 


(6.4) 


The total helical velocity U is, of course, a function of the radius. 

It is convenient to make the following assignments which are 

1 2 3 

consistent with the notation in Appendix A: x = x, x = y, x = z for 

12 3 

the Cartesian coordinates and y = y, y = r, y = £ for the helical 


coordinates. Using these replacements the next section will discuss the 
proper choices for A(r) and B(£). 
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r y 2 -CURVE 
j r y 3 -CURVE 



FIGURE 6.2. - PERIODICAL HELICAL COORDINATES 
FOR A CASCADE SHOWN WITH CARTESIAN COOR- 
DINATES. 
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6.3 Tailoring the Periodic Helical Coordinates 

The functions A(r) and B(£) are used to tailor the helical 
coordinate system to two sets of the streamwise helical sheets. The two 
sets of sheets of constant £ are illustrated in Fig. 6.3. They are 
evenly spaced in the circumferential direction. The first set of these 
sheets will contain the mean position of uncambered, symmetric, twisted 
blades. These correspond to the advance helicoids described earlier in 
Chapter II. The second set of helical surfaces is similar and is 
chosen such that each sheet lies midway between two neighbors of the 
first set, thus forming an alternating arrangement of periodically 
spaced helicoidal surfaces. The boundary conditions for the airfoil 
surfaces are applied very near members of the first set, whereas, the 
periodic conditions are enforced on the second. It should be mentioned 
that for asymmetrical blades or blades with camber, the mean blade 
position will not quite coincide with a helical sheet. The mean blade 
positions are assumed to lie near the first set of helical sheets so 
that small disturbance boundary conditions can be accurately applied. 

As mentioned, the helical sheets are surfaces of constant £. It 
is desirable for satisfactory application of blade surface boundary 
conditions that the ^-coordinate be orthogonal to the mean blade 
surface. In addition, simple handling of periodicity requires that 
there be no net change in axial distance when traversing a 
^-coordinate line from one periodic boundary to another. These two 
objectives can be met by properly choosing the functions A(r> and 
B(£) . 

Recall that the functions A(r) and B(£) enter the transforma- 
tion through the relationship given in Eq. (6.3), which is now written as 
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MIDCHANNEL MEAN BLADE 

HELICAL SHEET — v \ HELICAL SHEET 



FIGURE 6 . 3 . - ALTERNATING ARRANGEMENT OF HELI- 
COIDAL SHEETS FOR AN EIGHT-BLADED PROPELLER. 
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z = -y g + A(r)B(£) • 


(6.5) 


This specifies the axial coordinate in terms of the helical variables. 

It will be shown that by choosing suitable forms for A(r) and B(£), 
two of the metric tensor components anc * 023 vanish concurrently 

at special values of £. The y and i curves will be orthogonal 
where g ]3 = 0 and, likewise, the r and £ curves will be orthogonal 
where g 23 = 0. When both are zero, the E curves will be normal to 
the surfaces of constant E- The first objective is to arrange this 
to occur for those helical sheets containing the mean positions of the 
propeller blades. To assure that these surfaces are periodically spaced 
in the circumferential direction, B(£) is expressed as the product of 
two functions of E in the form 


b<c> = e<£) — 

tt 



( 6 . 6 ) 


where E is a damping function and E m is a positive constant 
establishing the period and is exactly half the distance between the 
blades . 

Meeting the conditions of orthogonality requires the inspection of 
the functional dependence of g^ 3 and 023' This will reveal what the 
final forms of A(r) and B(E> must be. First, the metric tensor 
component g^ 3 will be discussed, and then g 23 where additional 
constraints will be imposed on the axial velocity V(r). 

From Appendix A Eq. (A28), we have 

<6 - 7 


where the prime indicates differentiation with respect to the indicated 
argument. By using Eq. (6.6), this becomes 
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9 1 3 “ U R U Atr) 
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”~ s,n (l”)* E<5> cos (H] 


(6.8) 


As mentioned, E(E> is a damping function and is needed to assure 
a valid coordinate transformation. Useful selections of E(E> include 


E = exp(-a 




1I1\ 

*m ) 


, and E = exp|-a 



as well as polynomial forms. 


The significance of the bar above E in the exponential arguments 
will be discussed later. For simplicity only the first case, namely 
E(£) a exp^-a Q . will be presented since the other cases follow 


s i mi larly. 

The derivative E' (£) i s 


E * <C> 


= ± f exp/-a 

% V s m / 


(6.9) 


where the minus (plus) sign relates to positive (negative) values of E- 
By substituting the expressions for E(£) and E'(£) into 
Eq. (6.8) , we obtain 


'13 


Or 

U 


- q A(r)exp 


u Jii\r C os/ 

'M* - stn/f-Ol 

v ° w [ \ 

) * V 5 ™ JJ 


( 6 . 10 ) 


with the same meaning as above attached to the sign notation appearing 
before the last term. 

The periodic significance of E m is now determined by introducing 
a positive integer N g and allowing the value of E/E m to range from 
-1/2N to 1/2N C . Here, N R signifies the number of blades in the 
cascade and may be either even or odd. Since the value E m is half 
the distance between neighboring blades, the blades will be located at 
intervals of twice E . The total distance from - 1 /2Ng to l/2Ng 
equals one traverse around the cascade. By arbitrarily specifying 
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one blade to be positioned at E = 0, the blade locations take on the 

following values of £/£ . 

^m 

5 b 0, ±2, ±4, .... ±(N b - 1); for N b odd 

IjjJ = 0, ±2, ±4, .... ±<N b - 2), N b ; for N g even (6,11 

where the subscript B denotes a blade location. Also, at this point, 
the meaning of the over-bar above £ is made clear by defining J 
to be the value of E as measured from the nearest blade. Thus the 
extremes of E will be ±E m - With these definitions Eq. (6.10) 
for the metric tensor component can be written 



for values at the blade stations. However, since Eg = 0 from the 
definition of E and since cos^^ irj = 1 and sin^^ trj = 0, the 
value of the component g^ at the N g blade stations reduces to 


'13 




_ Or r V * / p \ 
'UR U 


(6.13) 


B 


Therefore, g 13 can be forced to equal zero at these periodic positions 
by simply prescribing that 

A<r) -r(s) 2 - <6 -' 4> 

The expression for the axial coordinate in Eq. (6.5) can now be 


expl ici tly written as 


V 


^m OR r 2 


2 - -y u + r r ^ exp 




(6.15) 


As mentioned above, the exponential factor serves as a damping function 
to guarantee that the coordinate transformation is well behaved; i.e.. 
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the coordinates should not fold nor should the Jacobian of the 
transformation vanish anywhere within the region of interest. In this 
regard, the quantity a Q is a nonnegative constant which increases 
in magnitude as the ratio V/QR = X decreases. For a given X, the 
larger the damping the more the ^-coordinate is forced to assume a 
purely circumferential direction. 

The expression for the metric tensor component g£3 can be 
written using Eq. (A30) in Appendix A as 
flr r U ' ( r )y 


’23 


U R U 


^ U ' (r) - + A'(r)B<£> A(r)B'(p. 


(6.16) 


From Eq. (6.6) and the values of the term A'(r)B(£) = 0 at the 
Ng blade locations. Furthermore, since g^ is zero at these 
positions, inspection of Eq. (6.7) reveals that AB 1 = (QR/V)(r/R) at 
the blades also. Then, for this case, the metric tensor component g 23 
reduces to 


g 23 




fir r U'y 
s - U R U 




(6.17) 


or m 2 ih 

= ’ v UJ u • 

By inspection, requiring V to be a constant is sufficient to render 
g 23 = 0. In this way, both g £3 and g ]3 are made equal to zero 
at the blade stations defined in Eq. (6.11). 

Thus, the E-coordinate lines will be orthogonal to the helical 
sheets containing the mean camber locations of the blades provided that 
A(r) and B(E) are defined as above and that the axial velocity V is 
a constant. The last constraint is a realistic restriction because a 
constant value of V is the most reasonable case and the one of most 


1 nterest . 
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The second objective is the requirement that the axial distance be 
the same for corresponding points at each periodic boundary. This is 
easily verified as being satisfied by inspecting the transformation of 
z on the helical sheets lying midway between successive blades. From 
Eq. (6.15), it is seen that 


Z PB " z 


V 

«PB = " Y ° 


(6.18) 


when 


£ ^PB 


’m 


’m 


±1, ±3, .... ±(N 
±1, ±3, .... ±(N, 


B 


2) , N b ; 

3) , ±(N, 


for Ng odd 


1 ) for 


Ng even 


B ’ _ B 

where the subscript PB denotes "periodic boundary." This relation 
shows that the net axial distance is not changed upon a complete traverse 
via a ^-coordinate line from one periodic sheet to any other. In 
fact, moving a value of twice £ m along any ^-coordinate restores the 
original axial location. That this is true can be seen by inspecting 
Eq. (6.15). The consequence of this is that any set of sheets separated 
by a value equal to the blade spacing 2£ m could be used as periodic 


boundary sheets. 


6.4 Final Form of the Periodic Helical Coordinates 

For convenience, the transformation of the helical coordinates 

12 3 12 

Y =y,r=y,S=y to the Cartesian coordinates x = x , y = x , 

3 

z = x is given in terms of the nonsuperscr i pted variables: 


x = 


y = 


r Sln (u Y + r) 
r cos (u Y + r) 


(6.19) 

( 6 . 20 ) 


2 ~ Y U + ir V 


^m QR r 2 


!xp (- a ° f) sin (H 


( 6 . 21 ) 


This is the set of coordinates which is used in the solution of the 
small disturbance equation presented in Chapter IV. Various other cases 
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of these coordinates could have been specified. For example, if either 
A(r ) or E(E> of equation (6.3) is chosen equal to zero, then the 
E-curves become circles about the z-axis. As another special case, 
if B<E) = E is chosen, then the E-curves are helices. On each 
helix the coordinate value varies as the arc length divided by the 
value of the radius at which that given helix lies. This distinguishes 
them from the original helical coordinates given in Appendix B. The 
metric tensors are given in the latter part of Appendix A for these 
special cases. 


VII. NUMERICAL APPROACH 


This chapter presents a general description of the numerical 

approach used to solve the small disturbance equation derived in 

Chapter IV. For the sake of convenience, the small disturbance equation 

1 2 3 

will be rewritten letting y = y, y = r, and y = as 
discussed in Chapter VI. With this notation Eq. (4.39) can be written 



3 <p _ <L 

" i <~\ r.1. — <>. . 


1. A. A, 


4 2 
3r 


5 3^ 2 


. a <L_¥_ . A + A 9 T 

6 3y3r 7 3y3E 8 3r3£ 


(7.1) 


where 


2MM-|-e 
A 1 = 72 / 3 “ 


(7.2a) 


A 2 - - | (H + 1 )MM J g 1 1 


(7.2b) 


J1 M 2 

A _ 9 

3 " .2/3 


(7.2c) 


e 22 
.2/3 9 


.2/3 


(7. 2d) 


a 5 = g 


33 


A - 12 

A 6 " _2/3 9 


(7.2e) 

(7.2f) 


2 1 3 

a 7 = 7T73 9 


(7.2g) 


A - q 23 

A 8 " .1/3 9 ‘ 


(7.2h) 
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In transferring the two terms on the left-hand side of Eq. (4.39) to the 
right-hand side in Eq. (7.1), they have been lumped together with the 
understanding that the variation of g 11 with respect to y is small. 
Recall that M is a function of r only and, therefore, the only 
dependence of or A^ on y is through g^ . 

It is convenient to define the two bracketed terms in Eq. (7.1) 
above by F so that: 



Equation (7.1) is then written as 


a 2 
O <p 


F + A, £-?■ + A c + A c fdb + A 7 + A a fjf 
Y 


4 7T T n 5 ac 2 + n 6 3y3r + n l 3y3^ n 8 3r3£ 
dr di; 


(7.4) 


n l 3y3t 

Except for the presence of the last two terms, this equation has the 
same form as the equation solved in Isom (Ref. 5) where an ADI method 
based on the Douglas-Gunn algorithm was used to solve the finite 
difference form of the equation for flow about helicopter rotors. The 
additional cross-derivative terms will be handled by generalizing the 
Douglas-Gunn algorithm. With modest changes this allows the numerical 
code developed for the helicopter problem to be used for the current 
work. 


It should be pointed out that Eq. (7.4) remains expressed in terms 
of unstretched physical variables. No mapping of the coordinate system 
has been carried out so as to produce a nonuniform grid. For 
simplicity, the numerical algorithm will be presented for the case where 
there is no coordinate stretching. Following this, the method of 
introducing coordinate stretching will be explained. 


Ill 


7.1 ADI Douglas-Gunn Algorithm 

For three dimensions, the ADI technique involves splitting the 

given equation into three separate finite difference equations which can 

be solved successively to complete one time-step increment. A current 

n n 1 

estimate of 9, say, 9" is advanced to <p n+ through two 

* * * 

intermediate values, which will be denoted by 9 and 9 , to complete 

a single stage of iteration. 

n * 

To begin an iteration at time step n, 9 is advanced to 9 

by solving the first equation with only the y-direction being 

differenced implicitly; this is called the y-sweep. Next, in the 
* ★ * 

r-sweep, 9 is advanced to 9 by solving the second equation with 

* * 

only the r-direction being implicit. Finally, in the £-sweep, 9 is 
advanced to 9 n+ ^ by solving the third equation with the ^-direction 
being the only implicit direction. 

The three equations for the respective sweeps are given for a 
uniform grid as: 

Y-sweep: 


At S y (<p ' * > = D y F I + Vrr* 


+ A 5 6 ^ 9 


+ Vyr* 


+ AyS^9 n + ( 7 . 5 a) 


r-sweep: 

Ai ** n ^4 ** n n 

— S(cp - 9 ) = D P T + 7- 5 (cp +<p) + ArSrt-9 
At y T T yI2rr T T 5 

/\ 

6 f ** n, A c n A f n 
+ 2 ~ V <Cf + 9 } + A 7 S y^ + A 8 S r^ • 


( 7 . 5 b) 
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S-sweep: 
A 

At v y 
A 


JT s „ (< f >n+1 - = D /l + Y 6 r r <9 + < p n> + Y S ?j: <c P n+1 + < f ,n> 


6 ** 
2 ~yv 


2Z s ( n+1 

2 V 


+ ^ 6 (<p + 9 ") + ^ S^(?" T ' + ?") + T* S. c <9 + 9 >• 


8 £ , n+1 

2“ S rC 


(7.5c) 


where is a special difference operator to be explained below and 
S rr’ S r£ are standarcl difference operators in the indicated 

directions. Here, F represents the bracketed terms of Eq. (7.1) which 
correspond to the flux in the Y - direction, and D^Fj 9^ ves the finite 

i_ L. 

difference approximation to 3F/3 y at the I-node for each set of 
values of J and K. The nodal values I, J, and K are associated with 
the y. r, and £ directions, respectively. The three directions have 
uniform step-sizes Ay , Ar, and A£, which may be distinct. 

An example of how the difference operators in Eq. (7.5) are defined 
in terms of difference approximations is given as follows 
n 1 


6 <p 

rr T 


(Ar )‘ 


r n - n n \ 

K.j+i.k " 2<p i,j,k + 9 I,J-1,kJ 


(7.6) 


An analogous expression may be readily written for S^. In the case of 
the mixed second-order operators, the following is used. 


r n J 

Yr 9 AYAr 


r n n n n \ 

vl , J+l ,K " 9 I f J,K " 9 I-1,J + 1,K + ^I-IJ.kJ • (7 ’ 


7) 


Corresponding expressions hold for 6^ and & r ^. 

The nonlinear term contained in D Fj is linearized by averaging 
at the n and * time levels by defining 

F I + 1 / 2 = A 2(Vl.J.O(Vl.J.O + * A 3 S y( 9 I,J,K + 9 I,J,k) <7 ‘ 

to be the flux at the midpoint of the I th -cel 1 on any node line given 
by J and K. The Murman-Cole type-dependent differencing scheme is 


8 ) 
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introduced to provide stable differencing as explained in Chapter III by 
defining 

D y F I = Ay [ e I (F I+l/2 " F I-l/2 ) + (1 e I-l )(F I-l/2 " F I-3/2 > ] * <7,9> 

where 

(1, V > 0 

e T = { (7.10) 

( °* V c < 0 

and 

V C = Vy + A 3‘ <7.11) 

This switches the difference equations at each grid point according to 
whether the flow field at that point is subsonic, sonic, supersonic, or 
a shock. 

A more convenient set of equations for numerical computations is 
obtained from the set given in Eq. (7.5) by subtracting Eq. (7.5a) from 
Eq. (7.5b) and Eq. (7.5b) from Eq. (7.5c), producing, with some 
rearrangement, the following set: 
y-sweep: 



+ A 7 S yE cpn + A 8 S r^ Cf>n; (7.12a) 

r-sweep: 

* n 
(9 - 9 ) ; 

(7.12b) 

£-sweep: 

*1 & _ ^5 ^7 — £ V n +l *\ 

At S y 2 2 6 yE 2 VeJ 9 “ 9 } 

A c A t A 0 \ ** 

r + r Sc + r s nJ (<p ~ 9 >• 



' A 1 A 4 

At \ ■ 2 + 6 rr 


- 2 ~ 6 yr) (<f ‘ 9 } = 


rr 




(7.12c) 
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The equations in Eq. (7.12) have been arranged into the so-called 
delta form where the unknowns on the left-hand sides are expressed as 
differences in the potential. This form has superior numerical 
properties compared with those in Eq. (7.5). The solution to the above 
set of equations involves no more than solving tridiagonal matrices, 
except in the case of a shock point where a quadridiagonal matrix occurs 
in the y-sweep. The potential at time level n+1 can be found from 
the potential at level n by adding the solutions for the delta 
differences from all sweeps to <p n as 

? n+1 = cp n + (cp* - cp n ) + (<p“ - 9 *) + W n+1 - ?**>• <7.13) 

The sequencing of the solution along lines of grid points for each 

of the three sweeps is illustrated in Fig. 7.1. The y-sweep is marched 

from upstream to downstream along rows of grid points where each row is 

characterized by constant values of r and E- The iteration in this 

sweep proceeds along rows in a plane of constant E by advancing 

sequentially to the row with the next higher value of r. When all rows 

of the current £-plane are completed, the process is repeated for the 

plane with the next higher value of E- When every plane has been 

* 

swept through, the solution has advanced to 9 . Similarly, the 
iteration proceeds in the r-sweep along rows of constant y and E by 

sequencing the rows through successive planes of constant E- This 

* * . , 
takes the solution to 9 . Since the value of E varies in the 

E-sweep, this last sweep proceeds along planes of constant r. When the 

last plane has been completed, the solution has advanced to 9 P+1 • This 

completes one full iteration. 

As is common in ADI methods, the time step is varied from iteration 
to iteration over the course of the calculation from some maximum value 
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At to some minimum value At. . This is done to improve the 
max min 

convergent rate of the calculation. The geometric sequence 

\(i-1 )/(N-1 ) 

1 = 1,2, 3, ...N (7.14) 

min/ 

is used for N iterations (N = 8) and then repeated until the total 
number of iterations has been reached. The total number of iterations 
is determined by a preset value for the maximum number of iterations, 
or by either satisfying a convergence criterion or exceeding an error 
bound. The range in the time step for the sequence of iterations 
addresses both high- and low-frequency components of the error. In 
general, the range of At mip and At max must be determined by 
trial and error and is strongly influenced by the size of the 
computational mesh. 

7.2 Grid Stretching 

Thus far, the physical mesh has been considered to be a uniform 
grid. However, it is preferable to have the grid points clustered in 
regions of high gradients and sparsely distributed in regions of low 
gradients. Grid stretching is a means of accomplishing this. It is 
used here to distribute the physical mesh points such that they are more 
heavily clustered near the airfoil than away from it with the greatest 
concentrations near the leading and trailing edges. The grid is 
smoothly stretched from the airfoil surfaces to a coarse grid at the 
outer boundaries of the flow field. The stretching is performed in all 
three coordinate directions. It is defined in a general sense as a 
mapping of the physical space yr£ to a uniform computational space 
Y r£ by 

y = y(y,r), r = r(r), £ = £(£)• (7.15) 


'At 

= ^mi n l At 
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It should be noted that the physical coordinates r and £ are 
stretched by their respective computational coordinates. However, y is 
a function of both y and r. The added dependence of y on r is 
necessary to accommodate swept wings. 

The coordinate stretching is introduced through the conventional 
chain rule formulas which are given in Appendix C. By replacing the 
various partial derivatives in Eq. (7.1) with those obtained by the 
chain rule, this equation is generalized to stretched grids. The use 
of coordinate stretching complicates, but does not change, the basic 
form of the ADI algorithm. 



VIII. DISCUSSION AND RESULTS 

This discussion focuses on the results obtained for flow over a 
single blade of an eight bladed cascade. As explained in the introduc- 
tion, the blades of this cascade are simple bicircular arc profiles with 
a maximum blade thickness of 5 percent. Furthermore, the planform of 
the blades is rectangular with the required spanwise twist made about 
the leading edge. This produces a blade of constant chord length. The 
aspect ratio of the blades is defined as the ratio of the blade-tip 
radius, as measured from the axis, divided by the chord length. All the 
results here are for blades with an aspect ratio of 4:1. The hub of the 
propeller system is placed at a radius of 0.375 R, where R is the tip 
radius; this gives an effective aspect ratio of 2.5 for a blade length 
measured from the hub, rather than from the axis. 

There are a number of reasons for using the blade geometry 
described above. The main reason is to provide a simple propeller 
geometry. By using a simple geometry, the flow will be uncomplicated 
by the effects which would otherwise arise by using a complex blade 
shape. The propeller geometry is simple for any specified advance 
ratio. Another reason is that the bicircular arc profile is widely used 
as a model profile in flow simulation studies. The front-to-back and 
top-to-bottom symmetry of the blade produces symmetrical flow for the 
case of a very high advance ratio and a zero angle of attack. For the 
case of a nonzero angle of attack, lift can be introduced simply by 
rotating the blade about the pitch change axis. Thus, a host of 
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propeller geometries can easily be generated by this method of 
specifying the blade construction. This, in turn, means that many 
different flow cases can be specified simply by prescribing the advance 
ratio, the angle of attack, and the Mach number of the approaching flow. 
Furthermore, any investigator who wishes to repeat the calculations can 
unambiguously reconstruct the propeller geometry. Because of these 
features, such a propeller system serves as a good prototype for flow 
Investigation. 

Four separate test cases are solved numerically in this 
investigation. These test cases are indicated in Table 8-1 by their 
respective values of the advance ratio X, the helical Mach number of 
the approaching flow M R , and the angle of attack a. The first case 
to be studied is the case of very high advance ratio X = 100 with a 
subsonic free-stream Mach number M R = 0.8 and a zero angle of attack 
a « 0.0°. Because this flow is essentially an axial flow, it represents 
a limiting case for zero propeller rotation. The second case studied is 
for an advance ratio X = 1, but with the same free-stream Mach number 
M„ = 0.8 and the same angle of attack a = 0.0° as used in the first 
case. The value X = 1 is typical for a propeller. The third case 
studied is for an advance ratio X = 1 and an angle of attack 
a = 0.0°, but for a transonic free-stream Mach number M R = 1.1. The 
values of X and M R for the third case are typical of an advanced 
turboprop. The fourth case is for an advance ratio X = 1 and a free- 
stream Mach number M R = 1.1, as in the third case, but now a nonzero 
angle of attack a = 2.0° is specified. Thus, the four cases present 
two values for each of the three varied parameters. The third case will 
include a separate study on the effect of grid refinement on the 
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solution. In addition, for the first three cases, the results obtained 
from the helical small disturbance (HSD) code are compared with results 
obtained from an Euler code (Ref. 6) These results are presented in the 
form of constant Mach number contours on the blade surfaces and selected 
cross-channel and blade-to-blade surfaces. All solution computations 
were made using a Cray X-MP/24. The number of numerical iterations used 
in computing the solutions for all test cases was 5000 for both codes. 
The computational time used by the HSD solution code was approximately 
60 percent of the computation time used by the Euler solution code. The 
Euler solution code required approximately three times the amount of 
computer memory of that used by the HSD solution code. Before 
discussing the results, the grids used in the computations for the HSD 
solution code and the Euler solution code will be detailed. 

Mesh lines for the HSD solution computations are shown (for the 
case X = 1) in Fig. 8.1 where, for clarity, only every third line is 
included from the leading edge to the trailing edge of the blade. A 
uniform grid extended over the blade surface and a stretched grid 
extended over the following regions: from the leading edge to the 

upstream boundary, from the blade tip to the outer radial boundary, and 
from the blade surface to the periodic boundary. The number of grid 
intervals for each computational direction is indicated in Table 8-2. 

The grid contained 30 intervals along the blade in the streamwise 
direction and an additional 11 grid intervals both upstream and 
downstream of the blade, for a total of 52 streamwise intervals. In the 
radial direction, 20 grid intervals were used, with half of these being 
on the blade. In the circumferential direction, 20 grid intervals were 
used from the lower to the upper periodic boundary. The total grid. 
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(A) VIEW NORMAL TO AXIS SHOWING PROJECTED 
VIEW OF HELICAL SHEET CONTAINING BLADE. 



(B) CROSS-CHANNEL VIEW AT LEADING EDGE. 



(C) BLADE-TO-BLADE VIEW 
ACROSS TWO CHANNELS; 

NOT ALL GRID POINTS ARE 
SHOWN ALONG BLADE SURFACE. 


FIGURE 8.1. - HELICAL COORDINATES USED IN SMALL DISTURBANCE COMPUTATION. 


122 


therefore was 52 by 20 by 20 intervals. The upstream and downstream 
boundaries were- located 4 chord lengths from the leading and trailing 
edges of the blade, respectively. The outer radial boundary was placed 
one blade diameter from the axis. 

This grid was used as the standard grid in computing most HSD 
solutions for the test cases. As discussed earlier for case 3, a series 
of grids in addition to this standard grid was used to Investigate the 
effect of grid refinement. This additional series of grids is described 
and the effect of grid density on the solution is discussed following 
the discussion of the case 3 results for the standard grid. 

For the Euler equation computations, the same unstretched grid 
extended over the blade surface. Otherwise, the grid was stretched as 
above, but with a different stretching function and no stretching in the 
circumferential direction. A representative grid, coarser than what was 
actually used for the Euler calculations, is shown in Fig. 8.2 for the 
case of X = 1. The number of grid intervals for each direction is 
given in Table 8-2. The position of the blade is indicated by the 
narrow opening (visible in Fig. 8.2(b)) of the grid lines near the hub. 
The grid of Fig. 8.2(c) illustrates how the streamwise grid transitions 
from an axial direction upstream of the blade to a helical direction 
near the blade and back to an axial direction downstream. Since the 
blade-to-blade direction is purely circumferential, this results in a 
high degree of coordinate shearing at axial locations near the blade. 

In addition, the chordwise distribution of grid points is not symmetric 
front to back along the blade surface, nor is it symmetric from the 
suction to the pressure side; this asymmetry increases with blade 
thickness and stagger. The upstream and downstream boundaries were 
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(A) VIEW NORMAL TO AXIS. (B) VIEW ALONG AXIS SHOWING THE MIDCHORD (C) VIEW Of BLADE-TO-BLADE SURFACE 

STATION. SHOWING THREE BLADES. GRID SHOWN 

IS COARSER IN AXIAL DIRECTION THAN 
THAT USED. 

FIGURE 8.2. - COORDINATES USED IN EULER COMPUTATION. 
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placed at two chord lengths from the leading and trailing edges of the 
blade, respectively. These boundaries were positioned differently from 
the corresponding HSD boundaries because of the difference between the 
two meshes in the streamwise direction. The outer radial boundary was 
located as in the HSD mesh, at one blade diameter from the axis. 

The computational results for the four propeller cases outlined 
above will now be discussed. Although the grids used in computing the 
solutions for these cases are coarse, based on two-dimensional 
standards, they are realistic for three-dimensional computations. The 
grid densities used in these calculations are sufficient to provide good 
overall prediction of the physics of the flow, but with some lack of 
detail, such as smearing of a shock. It is hoped that the following 
cases can serve as test cases for other researchers. 

Case 1) X = 100, M R = 0.8, a = 0° 

This case was chosen to examine the effect that blade cascading has 
on the solution. For the value of X = 100, the flow is essentially 
axial. Since the blade is symmetric from front to back and from top to 
bottom and a = 0, the solution should reflect this symmetry if no 
losses occur in the flow field. For the value of M R = 0.8 and the 
thin 5 percent thick blade, no shocks were detected in the flow field. 
The expected symmetry is noticeable in the solutions of both the HSD 
and the Euler codes. Figure 8.3 shows HSD Mach contours on the blade 
surface with the minimum contour being 0.75 and the maximum 0.9; the 
results are identical for the pressure surface and the suction 
surface. The results reveal the expected drop-off in Mach number 
with increasing radius. Similar results are shown in Fig. 8.4 where 
identical contour values are plotted for the Euler code. The contour 
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TIP 


L.E. 


HUB 
T.£. 



(A) PRESSURE SURFACE. (B) SUCTION SURFACE. 

FIGURE 8.3. - MACH CONTOURS OF SMALL DISTURBANCE COMPUTATION 
ON BLADE SURFACES: ADVANCE RATIO = 100, M R = 0.8, Ct = 0°. 
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shapes obtained from the two codes are very similar, with the only 
essential difference being that the Mach number produced at a given 
blade location is higher for the HSD code. 

Symmetry is again found for Mach contours in the cross-section 
planes given in Figs. 8.5 and 8.6 for the HSD and Euler codes, 
respectively. Again, the shapes obtained by the two codes are similar, 
with the HSD results showing more flow acceleration through the passage. 
The similarity in shapes indicates that qualitatively the solutions are 
being calculated correctly within the interior region of the flow, as 
well as at the blade surfaces. 

The solutions on the blade-to-blade surfaces are given in Figs. 8.7 
and 8.8 for the HSD and Euler codes respectively. In each case Mach 
contours are shown for three different span stations along the blade. 

The minimum contour value is 0.8. The values of the maximum contours 
are as follows: 1.0 for Fig. 8.7(a); 0.92 for Fig. 8.7(b); 0.86 for 

Fig. 8.7(c); 0.88 for Fig. 8.8(a); and 0.86 for Figs. 8.8(b) and 8.8(c). 
The results are symmetric and support the result that the HSD code 
predicts flows that agree with the Euler code except in magnitude, at 
least for subsonic and axial flow. 

Case 2) X - 1 , M R - 0.8, <x = 0° 

This case is presented to isolate the effect of blade rotation. 

The free-stream axial Mach number is only 0.5657, although M R = 0.8. 

The effect of operating at a low advance ratio is seen in Figs. 8.9 
and 8.10, which give Mach contours on the blade surface, for the HSD and 
Euler computations, respectively. These contours are given for constant 
Mach numbers that range from 0.6 to 0.8 in each case. Other than the 
expected result that the flow Mach number would increase toward the 



1-0.810 

2-0.826 

3- 0.843 

4- 0.860 

5- 0.876 

6- 0.893 

7- 0.910 


FIGURE 8.5. - MACH CONTOURS OF SMALL DISTURBANCE COMPUTATION IN 
CROSS PLANE AT MIDCHORD AXIAL LOCATION: ADVANCE RATIO = 100, 
= n 8 
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HUB TIP 


FIGURE 8.6. - MACH CONTOURS OF EULER COMPUTATION IN CROSS PLANE 
AT MIDCHORD AXIAL LOCATION: ADVANCE RATIO = 100, Mr = 0.8, 

Q = 0°. 



(A) r = 0.375 R. 


(B) r = 0.625 R. 


(C) r = 0.875 R. 


FIGURE 8.7. - MACH CONTOURS OF SMALL DISTURBANCE COMPUTATION ON BLADE-TO-BLADE SURFACES AT VARIOUS SPAN LOCATIONS: AD 

VANCE RATIO = 100, Hr = 0.8, a = 0°. 
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(A) r = 0.375 R. 
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1 - 0.80 
2 - 0.81 

3 - 0.82 
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5 - 0.84 

6 - 0.85 

7 - 0.86 



1 - 0.80 
2 - 0.81 

3 - 0.82 

4 - 0.83 

5 - 0.84 

6 - 0.85 

7 - 0.86 


FIGURE 8.8. - MACH CONTOURS OF EULER COMPUTATION ON BLADE- TO- BLADE SURFACES AT VARIOUS SPAN LOCATIONS: ADVANCE 
RATIO 100. Mr 0.8. Q = 0°. 
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6 - 0.766 
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(A) PRESSURE SURFACE. (B) SUCTION SURFACE. 

FIGURE 8.9. - MACH CONTOURS OF SMALL DISTURBANCE COMPUTATION ON 
BLADE SURFACES: ADVANCE RATIO = 1, M R = 0.8, a = 0°. 
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1 - 0.600 

2 - 0.633 

3 - 0.666 

4 - 0.700 

5 - 0.733 

6 - 0.766 

7 - 0.800 



(A) PRESSURE SURFACE. 


(B) SUCTION SURFACE. 


FIGURE 8.10. - MACH CONTOURS OF EULER COMPUTATION ON BLADE SUR- 
FACES: ADVANCE RATIO = 1, M R = 0.8, G = 0°. 
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blade tip, the Euler contours are asymmetric with the pressure surface, 
being nearly the Inverted Image (left to right) as compared to the 
suction surface. This might be the result of blade stagger which would 
give an Inverted Image for a symmetric blade. For the case of 
Isentroplc flow, the Mach number on the pressure surface at a given 
chord location would be the same as that on the suction surface If Its 
location was measured from the opposite end of the blade. To give an 
example, the maximum Mach number might occur at 60 percent chord, for a 
given span station, on the pressure surface; It would then have to occur 
at 40 percent chord on the suction surface. The reason for the observed 
difference In the magnitude between the pressure and suction contours 
for the Euler case is not known, but may, in part, be due to the grid 
asymmetry. While the results of the HSD contours are shown to be 
symmetric, there is no reason that the maximum Mach number must be at 
midchord. 

The blade-to-blade contours for this case are shown in Figs. 8.11 
and 8.12. The HSD results are given for the following range of Mach 
contours : 0.6 to 0.66 for Fig. 8.11(a); 0.66 to 0.71 for Fig. 8.11(b), 
and 0.75 to 0.78 for Fig. 8.11(c). The Euler results are presented for 
the same respective range of Mach contours. The primary difference 
between the two sets of contours is that the Euler computed contours 
more closely resemble contours about isolated blades. In the case of 
the HSD contours, they tend to shift upstream on the pressure side and 
downstream on the suction side of the blade so as to gradually join 
together at midchannel. 
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Case 3) x « 1 , M R « 1 .1 , a - 0° 

This repeats the previous case, except that now the free-stream 
Mach number is increased so that it has a value of 0.7778 on the axis 
and a helical free-stream value of 1.1 at the blade tip. The Mach 
contours on the blade surface are given in Fig. 8.13 for the HSD 
computation. The contours are shifted toward the trailing edge on both 
the pressure and suction surfaces, which show identical contours. Near 
the tip and trailing edge a very weak shock may exist. In the case of 
the Euler computation, the rearward shift of peak Mach number is more 
pronounced. A weak shock probably exists on the suction surface where 
larger gradients than on the pressure side are indicated in Fig. 8.14. 

The blade-to-blade contours for this case are shown in Fig. 8.15 
for the HSD computation and in Fig. 8.16 for the Euler computation. 

For both sets of results, the contour Mach numbers range from: a) 0.82 

to 1.0; b) 0.87 to 1.06, and c) 0.99 to 1.16. It is not clear that any 
shock exists for the HSD computation. However, a weak shock is 
observable in Fig. 8.16(c) of the Euler computation; it originates near 
the trailing edge of the suction surface and extends outward to a 
position upstream of the neighboring blade. 

The HSD program was also used to obtain solutions for this test 
case for three additional grid densities. The grids differed from the 
standard grid in the number of grid intervals used in the three 
coordinate directions, whereas the locations of the upstream, 
downstream, and radial boundaries went unchanged. Also, the type of 
stretching was the same as for the standard grid. The number of grid 
Intervals for each computational direction is presented in Table 8-3 for 
each of the three HSD refinement grids. These refinement grids are 
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FIGURE 8.13. - MACH CONTOURS OF SMALL DISTURBANCE COMPUTATION 
ON BLADE SURFACES: ADVANCE RATIO = 1, M R = 1.1, Q = 0°. 
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FIGURE 8.1A. - MACH CONTOURS OF EULER COMPUTATION ON BLADE SUR- 
FACES: ADVANCE RATIO = 1, M R = 1.1, 0 = 0°. 
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labeled coarse, medium, and fine, corresponding to their respective grid 
densities. For each of these additional grids, 6 grid intervals 
stretched from both the upstream boundary to the leading edge and from 
the trailing edge to the downstream boundary. Also, for each of these 
grids, 4 grid intervals stretched from the blade tip to the outer radial 
boundary. The variations among the grids were in the grid density used 
on the blade surface and the grid density used from blade to blade as 
follows: (1) the coarse grid contained 10 intervals from the leading 

edge of the blade to the trailing edge of the blade, 3 intervals from 
the hub to the blade tip, and 10 intervals from blade to blade; (2) the 
medium grid contained 20 intervals from the leading edge of the blade 
to the trailing edge of the blade, 6 intervals from the hub to the blade 
tip, and 20 intervals from blade to blade, and (3) the fine grid 
contained 30 intervals from the leading edge of the blade to the 
trailing edge of the blade, 12 intervals from the hub to the blade tip, 
and 40 intervals from blade to blade. 

Some solution results obtained using the refinement grids are 
provided in Figs. 8.17(a) to (c) where constant Mach number contours on 
the pressure surface of the blade are presented for the coarse, medium, 
and fine grids, respectively. These results are compared with the 
corresponding results presented earlier in Fig. 8.13(a) for the same 
range of Mach number contour levels as were obtained using the standard 
grid. In general, the contour patterns obtained for each of the three 
refinement grids are similar to the contour pattern for the standard 
grid. The highest Mach number contour is absent from the results shown 
in Fig. 8.17(a) since this level exceeds the maximum value computed for 
the coarse grid density. There are two effects of grid refinement to 
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(A) COARSE. (B) MEDIUM. (C) FINE. 

FIGURE 8.17. - MACH CONTOURS OF SMALL DISTURBANCE COMPUTATION ON PRESSURE 
SURFACE OF BLADE USING REFINEMENT GRIDS: ADVANCE RATIO = 1 , M R = 1.1, 
a o°. 
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be noted: (1) an increase in grid density results in a greater range 

of Mach numbers being obtained, and (2) an increase in grid density 
produces a more nearly symmetric contour pattern. Further refinements 
of the grid were made by increasing the density of the standard grid in 
each of the three coordinate directions, but only one direction at a 
time. Contour plots are not presented for these additional grid 
refinement results, but it is noted that they agreed closely with the 
contours provided in Fig. 8.13(a) for the standard grid. Based on 
these additional results and the fact that the contours of the fine 
grid refinement solution, presented in Fig. 8.17(c), agree closely with 
the contours for the standard grid solution, the refinement 
study indicated that the standard grid density is sufficient for most 
comparison purposes. 

To provide information on the convergence properties and the 
computational requirements of the HSD code, the following data are 
Included for the three grid refinement cases: (1) the overall reduction 

in the average solution residual; (2) the amount of computational time, 
and (3) the size of allocated computer memory. After computing 5000 
iterations on each grid, the following was found: (1) the average 

residual decreased by 7 orders of magnitude using the coarse grid, 3.5 
orders using the medium grid, and 2 orders using the fine grid; (2) the 
computational time required was 43 sec using the coarse grid, 204 sec 
using the medium grid, and 3204 sec using the fine grid, and (3) the 
memory allocated was 175,000 words for the coarse grid solution, 203,000 
words for the medium grid solution, and 330,000 words for the fine grid 


sol ution . 
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Case 4) X = 1 , M R = 1.1, a = 2° 

As a final case, the HSD code was used to recalculate the previous 
case except that a spanwise uniform angle of attack of 2° was used. The 
Euler code was not used for this case. The Mach contour plots are given 
in Fig. 8.18 for the blade surfaces and in Fig. 8.19 for the blade-to- 
blade surfaces. The effect of imposing an angle of attack on the blades 
resulted in a difference between the pressure and suction contours in 
the expected direction, i.e., the fluid velocity is now higher on the 
suction side. The blade-to-blade contours reveal that the flow is 
accelerated to a greater degree on the suction side. Although a weak 
shock may exist on either surface, no shock is noted to extend into the 
fluid from either surface of the blade. 
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FIGURE 8.19. - MACH CONTOURS OF SMALL DISTURBANCE COMPUTATION ON BLADE -TO- BLADE SURFACES AT VARIOUS SPAN LOCATIONS: 
RATIO = 1, Mfl = 1.1, a = 2°. 
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TABLE 8-1. - OPERATING 
PARAMETERS FOR THE 
FOUR PROPELLER 
TEST CASES 


Case 

X 

«r 

a 

■I 

100 

0.8 

0.0 


1 

0.8 

0.0 

n 

1 

1.1 

0.0 

H 

1 

1.1 

2.0 


TABLE 8-2. - NUMBER OF GRID INTERVALS IN EACH MESH REGION FOR THE HSD-SOLUTION 
STANDARD GRID AND BOTH EULER SOLUTION GRIDS 



Streamwise or axial 

di recti on 

Radi al 

di recti on 

Ci rcumf erenti al 
di recti on 


Upstream 

boundary 

to 

1 eadi ng 
edge 

Leading 

edge 

to 

trai 1 i ng 
edge 

T rai 1 i ng 
edge 
to 

downstream 

boundary 


81 ade 
tip to 
outer 
radial 
boundary 

Blade 

to 

blade 

Standard HSD- 
sol uti on 
gri d 

11 

30 

11 

10 

10 

20 

Euler solution 
arid for 
A = 100 

14 

30 

14 

10 

10 

20 

Euler solution 
cjrid for 

15 

30 

16 

10 

10 

20 


TABLE 8-3. - NUMBER OF GRID INTERVALS IN EACH MESH REGION FOR THE HSD-SOLUTION 

REFINEMENT GRIDS 



Streamwise direction 

Radi al 

di recti on 

Ci rcumferenti al 
di recti on 


Upstream 

boundary 

to 

1 eadi ng 
edge 

Leading 

edge 

to 

trai 1 i ng 
edge 

T rai 1 i ng 
edge 
to 

downstream 

boundary 


Blade 
tip to 
outer 
radial 
boundary 

Blade 

to 

blade 

Coarse HSD- 
sol uti on 
grid 

6 

10 

6 

3 

4 

10 

Medium HSD- 
sol uti on 
gri d 

6 

20 

6 

6 

4 

20 

Fine HSD- 
sol uti on 
grid 

6 

30 

6 

12 

4 

40 

































IX. CONCLUDING REMARKS 


In this thesis the flow over a propeller has been investigated. 

This investigation involved the following main elements: (1) using a 

potential formulation, the general tensor form of the equation governing 
the unsteady, inviscid, irrotational , and isentropic flow over a 
propeller in a noninertial, blade-fixed coordinate system was developed; 
(2) based on a coordinate system which is aligned with the undisturbed 
free-stream helical flow, a disturbance equation of equivalent accuracy 
to the general tensor equation was established in which the unknown is 
the perturbation potential, measured by its variation from the free- 
stream potential; <3) a systematic simplification of the perturbation 
equation was made, based on the scaling parameters characteristic of an 
advanced turboprop, thus leading to the establishment of a low- 
frequency, small disturbance equation for an unsteady, approximate (or 
small) perturbation potential; (4) the corresponding boundary conditions 
for the approximate perturbation potential were derived for both the 
solid surface and the farfield boundaries; (5) a new periodic helical 
coordinate system was introduced which provided for the straightforward 
treatment of the blade-surface boundary conditions and, also, the 
treatment of periodic boundary conditions for a cascade; (6) an ADI 
scheme, previously used to solve for flow about a helicopter rotor tip 
using Cartesian coordinates, was extended to include the capability of 
solving the propeller equation when more than one cross derivative term 
is retained; (7) a computer program, which was structured directly from 
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the program used to calculate the solution to the unsteady disturbance 
equation for flow over a helicopter rotor tip, was made which solved the 
unsteady small disturbance equation for propeller flow using helical 
coordinates; (8) steady-state solutions were calculated using the 
computer program developed here for four distinct propeller test cases, 
which included both subsonic and transonic free-stream helical flow, and 
(9) the results obtained for the four steady test cases were discussed 
and compared with corresponding results from an Euler solution program. 

The test cases presented in the last chapter were chosen because 
they serve as good prototypes for initial propeller studies. By 
inspection of the results obtained above, by either solution program, 
the flow of each case appears to lie within the range of computing 
capability of a small disturbance approach. In particular, no strong 
shocks were detected, nor did the computed flow depart by more than an 
acceptable extent from the free-stream state. Thus, for these studies, 
the flows appear to lie within the solution range governed, not only by 
the more complete Euler equations, but by the small disturbance equation 
as wel 1 . 

In regard to the above remarks, it is emphasized that an obvious 
goal in undertaking this investigation was to determine if a small 
disturbance approach, in the manner applied here, could be used to solve 
for the steady flow over lightly loaded propeller blades operating in 
the transonic regime. This goal was accomplished and illustrated for a 
number of test cases. Another central goal was to provide an estimate 
on the validity of the solutions obtained. Lacking exact solutions, an 
Euler solution code was used to provide comparison solutions for each of 
the test cases. The results of the test cases presented in Chapter VIII 
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indicate that, although both the helical small disturbance solution 
program and the Euler solution program provide reasonable solutions, 
enough differences exist between the two codes to warrant further 
investigation. These additional investigations should be conducted 
toward the purpose of verifying each of the programs, as neither can be 
considered as providing the exact solution. In particular, the 
difference occurring between the results of the two solution programs 
in the magnitude of the constant Mach number contours for case 1 should 
be resolved. Similarly, the difference between the results of the two 
programs for the blade-to-blade contours needs to be explained. Some 
of the differences may be attributed to the absence of a common mesh. 

As noted in Chapter VIII, the individual meshes used in these 
calculations were coarse. The density of the mesh was limited by the 
memory capacity of the computer. This limitation in capacity impacted 
the Euler program directly, since it required several times the memory 
capacity of the small disturbance program; the mesh density used in the 
small disturbance program was essentially set to match that used in the 
Euler program. Nevertheless, the grid refinement results obtained for 
the small disturbance solution of case 3 indicate that the grid density 
was reasonable, with the conclusion that further refinement of the grid 
was not going to result in large changes in the distribution of the Mach 
number contours. In addition to the mesh coarseness, other variations 
in mesh characteristics existed between the two codes. Some of the 
computational differences between the two codes might have been resolved 
if a common mesh had been used. As a final remark concerning solution 
verification, it is often the case that the verification of a solution 
is the most difficult step and one that is dependent on the existence 
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of solutions from other programs. In the case of the propeller, the 
capability to solve the flow is only now emerging, and hence the 
verification process will certainly be an ongoing one. It is hoped that 
the flow cases presented here will serve as test cases for other 
efforts . 

A few comments will now be made concerning the potential value of 
the small disturbance computer program. In its present form, the 
program is capable of solving flow about lightly loaded blades, where 
the flow is steady in a blade-fixed coordinate system. However, the 
governing equation developed and presented in Chapter IV is for an 
unsteady flow. The decision to study only steady flow cases in this 
investigation was made not because of a particular limitation in the 
method, but because the steady flow problem was deemed sufficiently 
important and difficult to be treated separately. The relative value 
of the small perturbation method, as compared to methods based on 
equations valid for more general flows, such as the full potential 
equation or Euler equations, is greater for unsteady flow than for 
steady flow. The reason for this is that the computational resources, 
both in the memory capacity and speed of the computer, required to solve 
an unsteady flow are much greater than the resources required for the 
solution of a steady flow. Clearly, the reduction in computational 
expense offered by a small disturbance approach in comparison to, say 
an Euler equation approach, will be compounded when the use of either 
computer storage or computer time is increased. 

As an example of a practical circumstance that illustrates the 
additional resources needed for solving an unsteady flow, the 
restrictions on the time step is an obvious case in point. The time 
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step must be restricted to a degree which allows the unsteady phenomena 
to be resolved; this is a significant limitation, and one which is often 
further complicated by the existence of fundamental time scales of 
widely differing magnitude. The minimum time step, corresponding to the 
smallest time scale required to resolve the flow, essentially determines 
the overall expenditure of computer time. 

As another example which illustrates the additional demands made 
on computer resources in the case of an unsteady problem, the 
computational mesh size is mentioned. For most steady problems 
involving a farfield boundary surrounding a body, the computational mesh 
is stretched, so as to produce a grid having larger intervals away from 
the body than near it. This is especially typical with respect to 
transonic flow problems, where the decay of perturbations traveling 
lateral to the flow is extremely small, and consequently, the farfield 
is often placed 100 body lengths away. With a distant placement of 
farfield boundaries, a stretched grid is highly desirable, as it offers 
a significant reduction in mesh points in comparison to an unstretched 
grid; the magnitude of the grid stretching, of course, affects the 
solution accuracy. However, when using a stretched grid for computing 
an unsteady flow, additional errors, beyond those encountered for steady 
flow, are introduced. An example of the cause of such errors is the 
inaccuracy that occurs in the replication of an outgoing wave when it 
travels through a grid which is stretched. Near the body where the grid 
is finely spaced, the characteristics of the wave may be well 
represented; away from the body, where the grid is coarse, the wave may 
be distorted, both in magnitude and frequency. Furthermore, the 
resulting distorted wave may be reflected back toward the body, 
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whereupon, after further distortion occurring on its return path, it 
later reaches the body, but with a different wavelength and a different 
magnitude than it would have had if a uniform grid had been used. 

Additional concerns, beyond the time-step and mesh size 
limitations, arise in the case of unsteady flow, and include the proper 
treatment of an aerodynamic body displaced from its equilibrium position 
due the action of fluid forces. Many of these concerns may be 
adequately addressed by adopting a perturbation approach. The example 
of the restriction on mesh size discussed above illustrates the 
advantages of a perturbation approach of the type used in the present 
investigation. For the same grid density, a small disturbance approach 
will require only a small fraction of the amount of computer memory as, 
say an Euler solution approach. With a given amount of computer memory, 
a finer grid, stretched or unstretched, can be used with the small 
disturbance approach in comparison with the same type of grid used for 
the Euler code. Furthermore, a small disturbance approach is computa- 
tionally more efficient with respect to computer time, in general, than 
methods solving the more general equations. Since, as was discussed 
above, the unsteady flow case requires a greater number of mesh points, 
the computational savings of an efficient code is compounded above the 
savings realized for steady flow calculations because the computational 
time required to complete a solution iteration increases in proportion 
to the total number of mesh points. In view of the above remarks, and 
considering current computer capabilities, the small disturbance 
approach presents a viable method of solving both steady and unsteady 
propel ler flow. 
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Based on this conclusion, a few suggestions are made regarding 
possible future investigations of propeller flow using the helical small 
disturbance approach. These suggestions include methods of improving 
on the above investigation and also possible extensions to the scope of 
the investigation. As a first step, test-case comparisons similar to 
those studied here should be computed using a common mesh. As a means 
of providing a common mesh, a more general helical coordinate system 
than the one presented in this thesis has already been developed. This 
newer helical coordinate system conforms to the exact blade shape, 
rather than the mean chord position. This means that the small 
disturbance boundary conditions can be replaced by actual blade surface 
boundary conditions. The new coordinate system is suitable for other 
inviscid flow solvers, and thus provides a common mesh upon which 
comparison solutions could be computed. Additional studies could also 
be made investigating the effect of using the exact blade boundary, 
rather than the mean chord position, as the location for the blade 
surface boundary condition. For unsteady flow which may involve blade 
flutter, the equilibrium position of the blade may be used in place of 
the exact location of the blade, if the blade deflection is small. In 
this way, the small disturbance approach to the treatment of the blade 
boundary condition can be accurately extended to unsteady flow. These 
are some of the possible areas which may be examined in future 
studies. Should any additional work along these lines be undertaken, 
the approach chosen should be as simple as reasonable accuracy allows. 



APPENDIX A - HELICAL COORDINATE TRANSFORMATION FOR A CASCADE 
In this section, the metric tensor components are derived and the 
Jacobian is determined for a general transformation between orthogonal 
Cartesian coordinates and helical coordinates y^ . The 

Cartesian coordinates will be given as functions of the helical 
coordinates which is expressed formally by 

x 1 = i.j = 1,2,3. (Al) 

Both x 1 and y 3 are right-handed triplets. To avoid carrying 
along the superscript notation, the following assignments are made 

x 1 = x, x 2 = y, x 3 = z (A2) 

for the Cartesian coordinates and 

y 1 = y, y 2 = r, y 3 = K (A3) 

for the helical coordinates. The transformation of coordinates is then 
defined through 

x = x(y,r,£) = r sin d <A4) 

y = y(y,r,£) = r cos i? (A5> 

z = z(y,r,£) = -y g + A(r)B(£) (A6) 


where 


0 + I 

" " U R 


(A7) 


and is measured from the y-axis as is shown in Fig. (6-2). The total 
or helical velocity U and the axial velocity V are each functions 
of r and are connected by the rotational velocity Q through the 
Pythagorean relation given by 
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U = U(r) = [(Or) 2 + V 2 ( r )J (A8) 

The significance of the functions A(r) and B(£) and their actual 
forms were considered In Chapter V. Here they are assumed to be 
reasonable functions of r and £ respectively. 

By defining 


u’ - fr • v ' - a? ’ A ' - a? * and B ' - f • <A9) 

the individual partial derivatives of the Cartesian coordinates, which 
will be needed, can be expressed as 


ax 

ay 


1 _ ax _ Or 

i - a Y " u cos 0 


(A10) 


9X ^ _ 3X , „ 9$ q 

gy 2 “ 9r " s ^ + 9r cos § 


(All) 


_ r 

g y 3 35 ‘ R 


(A12) 


E = - fp sin d 

ay 1 8 y U 


( A 1 3 ) 


^2 i §* . cos g - r If sin i 


(A14) 


ax _ 3y r , 

g y 3 = 35 - - R 


(A15) 


ax _ az 
ay 1 = ^ 


V 

u 


(A16) 


ax J , az 
ay 2 = 8r 


n~ )y + a 1 b 


(A17) 


3x_ _ az ad i 
3y 3 = * 


(A18) 
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A necessary and sufficient condition for the transformation 
relations to be independent is that the Jacobian J of the 
transformation, namely 

3x ] Sx 1 3x1 

1 2 3 

3y 3y 3y 

2 2 2 

3x 3x 3x 

1 2 3 

3y 3y 3y 

3 3 3 

3x J 3x 3x 

1 2 3 

3y 3y 3y 

be nonzero. The partial derivatives given by equations (A10) 
through (A18) can be used to evaluate J; which is expressed in mixed 
form to reduce the algebra, as 

cos 0 si n 0 + r cos 0 ^ cos i? 

J = - sin 0 cos 0 - r sin d - ^ sin 0 . (A20) 

- V 32 3z 

U 3r 3J- 

Expanding the determinant produces the following: 

J = jj£ || (cos 2 d - r cos d sin «?) + g £ (sin 2 d + r cos 0 sin 

+ ~ || (sin 2 t? + r cos sin 0) + g £ (cos 2 - r |p cos 0 sin tf) 

‘ S F fr cos ^ s1n 0 + £ (T fr cos d Sin * = 0 R + IT Ie ' <A21) 

From this, the transformation is seen to be valid in the domain where 



( A 1 9 ) 


J * 0, i .e. , 
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9z „ _ y_ r _ x 

di* QR “ A - 


(A22) 


The covariant metric tensor components can be determined using 


3x m ax m 

i a J ’ 


^ ay 1 ay 

This will be used to find each of the six distinct g. ^ . First, to 
find g^ the following three terms are easily developed: 


(A23) 


8x1 axl - / 3x x2 
ay 1 ay 1 

ax 2 ax 2 _ ( §y x2 
ay 1 ay 1 


Qr ^ cos 2 0 



axlaxl./azV /y\ 2 

v u / 


ay 3y 

Adding the three terms as prescribed by Eq. (A23) gives 

which follows from Eq. (A28). Similarly, the terms for g 22 are 


(A24) 


ax 1 ax 1 


ay ay 

2 2 
3x_ 3x_ 

2 2 

ay ay 



ax 


,_2 


dO 


= lffJ = sin 0 + 2r g cos 



3r 


3r 


-ii^u 


2 2 
3y^ 3y^ 


3r 



cos 2 0 


and, thus 


*22 


1 + r 


2 dd 



3r 



(A25) 
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For 9 3 3 


and, thus 
For g 12 


8x1 8x1 - / 9x\ 2 
ay 3 ay 3 = W 





■(r) 2 co;2 

" (r) 2 si " 2 J 

= (AB 1 ) 2 


g 33 ' 




ax 1 ax 1 _ ax ax nr . nr at? 2 , 

^2 = a^ a? = [T cos » sin * + u" r 9? cos * 


9x 2 8x 2 _ 9y <^y fir g . . fir ^ M tin 2 

8 1 8y 2 " 9y 9r " U C0US 'M + U r gr sin 


and, thus 


9x 3 9x 3 _ 9z 9z 


1 2 
9y 9y^ 


9y 9r 


[-5(v - a u ) - A B ] 


fir 9<? 


g 12 = U 9r ' " a - 


V 9z 
U 9r 


For g 


13 


9x \ 9x \ _ 9j< 9x 

ay' ay 3 = ^ 35 


fir r 2 a 

u“ R cos 0 


2 2 

9x 9x _ 9y 9y fir r .2 , 
ay 1 9y 3 = = U R 


9x 3 9x 3 - 9z 9z V flRt 
9y ] 9y 3 = 9 Y = " * 


(A26) 


(A27) 
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and, thus 


_ _ r V -p, 

9 1 3 — U R ' U AB ' 


(A28) 


For g 


23 


ax 1 3 x ] _ 3 x 3 x r , , , r dd 2 « 

^2 ^3 = 8? 8| = R C0S * Sin * + R f 3r C0S * 

3 x 2 3 x 2 _ 3^ 3y r , „ r 3i? .2 ^ 

^2 ^3 = 3? if = " R C0S * Sin ° R 3r S1 " 


3 x 3 3 x 3 _ 3z 3z 
3y 2 3y 3 = 8r 


(j2 U' -f)"*' 1 ]* 6 ' 


and, thus 


r 3z 3z 

*23 " R 3r + 3r 


or 


Or r U ' y 

g 23 “ “ U R U + 


(j2 u ' - r)v - A ' B ] AB> 


(A29) 


(A30) 


We can collect g^, g 22> 933 . 9 12 > 9] 3 * and 9 23 ^ nt0 the m ^trix 
1 


'Ij 


Or 8 j? V 3z 

U 8 r U 9r 


Or r V 3z 

UR U 8 £ 


Or 3i? V 3z 

U 3r ' U ar 


flr r V 3z 

L U R 0 3? 



r dd 3z 3z 
R r 3r + 3r 9^ 

rr\ 2 /3z \ 2 

CrJ + lad _ 


(A31) 


The Det g... can be readily calculated. To simplify the notation, 


let 


a = 


3z 

3r ’ 


6 - If • and n * r if 


(A32) 


where a, 13, and n are used only for convenience and are independent 
of their use elsewhere in the text. The determinant of g^, when 
expanded, gives 
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Further expansion and use of the identity 



yields 


(A34) 



Cancellation reduces Eq. (A35) to the following form 
or 

a n of n l x v 3zV 

9 - Det 9 ij = \R U + 0" M) • 


(A37) 
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Comparison of this result with Eq. (A21) shows that 

3 = <g>' /2 , 

exactly as It must. 

The contravari ant metric tensors g ig can be calculated 
and g using the general relation 

jj _s!l 


g ’ 


where are the corresponding minors of the determinant < 

g^ are determined as follows: 


( ! 


11 1 


g11 _ g l g 22 g 33 g 23 


i 2 2 

I + n + a 


) 


(ID * - R n + “®] | 

1 [(I) 2 . ♦ (jjfc - 2 { na B ] 


therefore , 


11 1 


(r) + ° 2 * (r a - nl3 )] 

) 

(lr r ■ o B ) 2 j 



9 1 1 g 33 g l 3 , 



therefore, 


g33 ■ g (^1 1^22 - »?j) 


1 + 


2 2 / Qr V \ 2 

n + a n - U a ) 


(A38) 
from g. j 

(A39) 
I. The 


(A40) 


(A41 ) 
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therefore, 



(A42) 


--ij(o E "-o«)[(ii) 2 *» 2 ]-(8 t i-S»)(iin + -)j 
■ ' g | [s 0 * T 5 ] n0 ' I [ij 0 * r °] “I 

therefore , 

9 ' 2 ■ yi(ii a - nB ) • <A43> 

9 = g (g 21 g 32 " g 31 g 22 ) 

9,3 ^[(pl“)(i"*' ll )-(n-0 C )0 * I 2 * a 2 )]- <A44> 

Finally 

9 = " g <g ll g 32 " g 31 g 12 ) 

■■j|(c)^*ros**(rs-(r) ! «]j 

1 (r V Qr \( V fir \ 

= - g 0 + cr 7lo n + 0" a ) 

therefore , 

23 1 / V fir \ 

9 ■ ‘ y3^u n * D“ a ) 


(A45) 
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The g^ can be expressed in the symmetric matrix form similar to 


Eg. (A 31 ) as 



) fe ) 2 * 02 T* (S * - n 0 ) 

9 [(r n - 0 a )(s n * aB ) 


* G 0 - nB ) ] 

- (r s- s«)o - 2 * D 2 )] 

g’ J - 

g' 2 1 

1 (\l Qr \ 

- yf Vu 11 + 0- a ) 

_ — 1 


g' 3 g 23 

g [’ *(o n * r a ) J 



(A 46 ) 


Expanding the determinant of the above matrix for the contravariant 
metric tensor components provides a check of the algebra and gives 

° et 9 ' J ' 7 | [® * 132 * (s “ * nl3 ) ][' * (0 n + r a ) ] 



This simplifies to the following form after cancellation of some of the 
terms 


Det g ij 



- [( n0 - s a )(o n + r a ) 

0 a )(s 11 * aB ) * (r s - 0 B ) n + n2 * a2> ] | • 



166 


In turn, this expression can be expanded to give 


™ - hW - (i ^ 


2 Clr „ r V Qr r 2 Or r 2 
+ U" naB " R 0 na " 0“ R a “ 11 ° n 


U R 


fir „ r V V 2 n Or r V n Or r 2 V 2 R 

\T naB + RD na + D aB + u r ~ u 13 u R n u nB 


ar r 2 V 2 R V1 
Ti* - u a a ) ] 


More cancellation brings 


Det 9 ,J - Jj [(D 2 * » 2 - (fj (O' * 


Qr r V 


U R U 


0 - 0 



Using the identity given in Eq . (A34), leads to 

v 2 


Det g 

Using Eq. (A37), we finally have 


j j . ]_ [^r yV 

' ‘ g 2 LV R U / 


0 Or r V n 
+ 2 U R U 13 + 
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Special Cases of Coordinates 

Two special cases are now given. In the first case, for A = 0 in 
Eq. (A6) and with V equal to a constant, we have the case where £ 
is a purely circumferential direction and 


1 



g 


13 



(A48) 
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In this case J = (r/RXV/U). 

In the second case, again, Visa constant, here, however, 
2 

A = (QR/V)(r/R) and B = This produces a coordinate system 
where the £ curves are helices with the arc length inversely 
proportional to r. The metric tensor components are 



and. 


(A49) 


(A50) 



In this case J = (r/R)(U/V) 



APPENDIX B - HELICAL COORDINATE TRANSFORMATION FOR AN ISOLATED BLADE 
The transformation is given here for the helical coordinates used 
for flow about an isolated blade with a tip radius R which is rotating 
with constant angular speed Q and advancing at a constant axial 
speed V. The covariant metric tensor components g^ and the 
contravar i ant components g^ 3 are also given for reference. Let 
xV.x 3 be a right-handed orthogonal Cartesian coordinate system and 
y\y 2 .y 3 be a right handed helical coordinate system as shown in 
Fig. (6.1). The transformation of the helical coordinates to the 
Cartesian coordinates is given by 


x 


1 



x 


3 



sin $ 
cos d 


♦ § y V 


(Bl) 

(B2) 

(B3) 


2 , 

where the angle $ is measured from the y -axis as shown in 
Fig. (6.1) and given by 

, 0 1 V ^ 

* U y ♦ u 2 • 

Here, the total velocity U is given as 

r 2 ^ 2*p ^ 

U = (Qy^> + V* 

The covariant metric tensor components g.^ are 


(B4) 


(B5) 


to 
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and the contravariant components are 

1 + ( 9i2 )2 _g 12 

g 1J = -g 12 1 

9 12 9 23 _9 23 

Since for these coordinates Det g. . E g = 
of the transformation is unity, J = 1. 


g 12 g 23 

" g 23 

1 + (g 23 > • 


(B7) 


then the Jacobian J 



APPENDIX C - PARTIAL DERIVATIVES FOR COORDINATE STRETCHING 
Given coordinates y, r, % which are functions of the coordinates 


y,r,E defined by 



Y = y(y, r), r = r(r), - 

(Cl) 

the following 

relationships between the partial derivatives are 


obtained: 

3<p 3<p 3y 

3y _ 9y 9y 

(C2) 


3cp 3<p 3y 3<P 3r 

3r = 3y 3r + 9r 3r 

(C3) 


3<p 9<g 3E 

3E ' 9C 9^ 

(C4) 


3 Z ip / §i\ 2 . <b£. 

a Y z ’ 3, Z W * 3y 2 

(C5) 

? 2 
9 <p 9 <p 

3r 2 3y 2 

(liY , 2 IS 3f , |£ £i „ 8% (|rV + H ai| 

\3r/ 2 3y9r 3r 3r 3y 8r 2 g-2 \3r/ 3r 9f 2 

(C6) 
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(C7) 


o 9 2- 
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